Calibrate ECMS Block for Objective Hybrid Vehicle Fuel Economy Assessment

Open Live Script

When you compare fuel economy changes in a hybrid vehicle drive cycle, you must account for the amount of energy stored in the vehicle battery. Otherwise, net battery charging or depletion during the drive cycle can artificially inflate or deflate the fuel economy estimate for the drive cycle. Powertrain Blockset™ allows you to simulate the fuel economy, performance, and emissions of hybrid electric vehicles with standard hybrid powertrain architectures P0, P1, P2, P3, and P4.

Powertrain Blockset uses an industry-standard equivalent consumption minimization strategy (ECMS) to optimize the fuel economy of hybrid reference applications by optimally trading off internal combustion engine (ICE) use versus electric motor use during a given drive cycle.

This script optimizes the weighting factor parameter, ECMS_s, automatically. The script runs a drive-cycle simulation, measures starting and ending SOC, and repeats the process until the difference between starting and ending SOC is minimized.

During the optimization process, an optimization plot indicates the progress of the optimizer toward net zero change in SOC over the given driven cycle for the chosen vehicle configuration.

The optimization may take more than one hour to complete. After optimization completes, the ECMS_s parameter for net zero change in cycle SOC is stored in the Optimal Control Model Reference model workspace.

The Equivalent Consumption Minimization Strategy block in the hybrid controller of each hybrid reference application optimizes trading off engine and motor use during the cycle. The Equivalent Consumption Minimization Strategey block uses a weighting factor parameter, ECMS_s, that you can set so that the hybrid battery state of charge (SOC) at the end of a given drive cycle is the same as the SOC at the beginning of the drive cycle. The block optimizes fuel economy.

Screenshot 2024-06-13 004522.png

Screenshot 2024-06-13 115828.png

ECMS Calibration Script

Use this code to calibrate the ECMS_s block parameter for net zero change in battery SOC for any of the Powertrain Blockset Hybrid Reference Applications P0, P1, P2, P3, and P4.

Choose HEV Reference Application

First, choose the HEV reference application. This example uses the P2 reference application.

% Open a given reference application
pxName="P2"; % HEV powertrain architecture configuration: P0,P1,P2,P3,P4
eval("autoblkHev"+pxName+"Start"); %Start a Px Reference Application
SOCinit = 0.6; % initial SOC, unit is in [0, 1]

Next, execute the optimization function Calibrate_ECMS_s.

As the optimization runs, you can view progress in the FTP75 Drive-Cycle SOC Balance plots. The dSOC versus ECMS_s plot shows the difference between the optimized SOC and the actual SOC for each weighting factor. As the optimization converges, the difference approaches zero. The SOC versus ECMS_s plot shows the actual SOC for each weighting factor.

ECMS_s=Calibrate_ECMS_s(pxName,SOCinit); %Run ECMS_s optimization and return optimal ECMS_s for net zero change in SOC over the drive-cycle

$Figure contains 2 axes objects and another object of type subplottext. Axes object 1 with ylabel dSOC contains 2 objects of type line. One or more of the lines displays its values using only markers Axes object 2 with xlabel ECMS\_s, ylabel SOC contains 2 objects of type line. One or more of the lines displays its values using only markers$

Elapsed time is 590.795579 seconds.
Search converged. ECMS_s parameter updated in model. 
ECMS_s = 4.442937, SOCEndTrg = 60.000000, SOC_end = 59.755828

ECMS_s Calibration Function

function ECMS_s=Calibrate_ECMS_s(pxName,SOCinit)

    % calibration_ecms_script
    battMdl = "BatteryMapped"; % battery model - for SOC sweep
    mdlName="ConfiguredHev"+pxName+"VirtualVehicle"; % Hev model
    ctrlMdl="Hev"+pxName+"OptimalController"; % Controller model
    SOCEndTrg = SOCinit*100;
    maxIter=4;
    
    % Set up optimization progress plot window size and position
    x0=600;
    y0=1040;
    width=600;
    height=280;
    
    tic
    load_system(mdlName);
    load_system(ctrlMdl);
    load_system(battMdl);
    blk = mdlName + "/" + "Scenarios"+"/"+"Reference Generator"+"/"+"Drive Cycle"+"/"+"Drive Cycle Source";
    m = get_param(blk, 'MaskObject');
    ECMS_CurrentCycle = m.Parameters(1,1).Value; % Current cycle setting.
    % get current ECMS_s value
    ddObj = Simulink.data.dictionary.open('VirtualVehicleTemplate.sldd\HevPxVCU.sldd');% find data dictionary
    dataSectionObj = getSection(ddObj, 'Design Data');
    ECMSName="CtrlVcuHev"+pxName+"ECMSs";
    ECMS_obj = getEntry(dataSectionObj, ECMSName);
    if  isnumeric(ECMS_obj); ECMS_CurrentValue = ECMS_obj; end
    if ~isnumeric(ECMS_obj); ECMS_CurrentValue = ECMS_obj.getValue; end
    
    % extract initial value/name of ECMS tuning parameter
    p = Simulink.Mask.get(ctrlMdl+"/ECMS");
    baseParamName=p.getParameter("ECMS_s").Value;
    ddObj1 = Simulink.data.dictionary.open('VirtualVehicleTemplate.sldd\BatteryDCDC.sldd');
    dataSectionObj1 = getSection(ddObj1, 'Design Data');
    battChrgMaxValue_obj = getEntry(dataSectionObj1, 'PlntBattChrgMax');
    if  isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj; end
    if ~isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj.getValue; end
    
    % set to a temp name
    set_param(ctrlMdl+"/ECMS",'ECMS_s',"ECMS_s_tune")
    save_system(mdlName,[],'SaveDirtyReferencedModels','on');
    save_system(ctrlMdl,[],'SaveDirtyReferencedModels','on');
    
    % Enable SOC recording
    x1_handles = get_param(mdlName+"/Visualization/Rate Transition1",'PortHandles');
    x1 = x1_handles.Outport(1);
    Simulink.sdi.markSignalForStreaming(x1,'on');
    
    % arrays used to store ECMS_s and SOC_end values
    xc = zeros (3,1);
    yc = zeros (3,1);
    
    % plot showing the ECMS_s and SOC_end
    subplot (1,2,1);
    set(gcf,'position',[x0,y0,width,height])
    xlabel({'ECMS\_s',' '})
    ylabel('dSOC')
    hold on;
    subplot (1,2,2);
    xlabel('ECMS\_s')
    ylabel('SOC')
    hold on;
    sgtitle([ECMS_CurrentCycle ' Drive-Cycle SOC Balance']);
    % set model workspace variables
    in = ModelSetVariable (mdlName, battChrgMaxValue, SOCinit,battMdl, SOCEndTrg, ctrlMdl);
    % get initial 3 ECMS_s points and SOC_end values
    [x, y, ECMS_s, SOC_end, solution_found] = dSOC_generate_starting_points (in, ECMS_CurrentValue, SOCEndTrg);
    
    if ~solution_found
        
        for i = 1 : maxIter
           
            % solve second order equation to get z = ECMS_s corresponds to y = SOCEndTrg

$\begin{array}{l} S O C E n d T r g = S O C_{i n i t} = a x^{2} + b x + c, where a, b, c satisfies y_{i} = a x_{i}^{2} + b x_{i} + c, 1 \leq i \leq 3 . with \\ x_{1} = {ECMS}_{s}^{(1)}, y_{1} = S O C_{e n d} ({ECMS}_{s}^{(1)}), x_{2} = {ECMS}_{s}^{(2)}, y_{2} = S O C_{e n d} ({ECMS}_{s}^{(2)}), x_{3} = {ECMS}_{s}^{(3)}, y_{3} = S O C_{e n d} ({ECMS}_{s}^{(3)}) . \end{array}$

            [z] = Second_order_roots (x, y, SOCEndTrg);
            in = in.setVariable("ECMS_s_tune", z);
            [SOC_end, dSOC] =dSOCsim_v1(in);
            subplot (1,2,1);
            plot(z,dSOC,'bs','MarkerSize',15)
            grid on
            hold on

            subplot (1,2,2);
            plot(z,SOC_end,'bs','MarkerSize',15)
            grid on
            hold on
            
            if (abs(SOC_end - SOCEndTrg) <= 1) % check if a solution is found
                ECMS_s = z;
                x(1) = ECMS_s;
                y(1) = SOC_end;
                solution_found = 1;
                break;
            end
            
            xc(1:3) = x (1:3); yc(1:3) = y (1:3); 
            x(1) = z; y(1) = SOC_end; % since z and SOC_end are new points, we use them and put into x(1), y(1)
            [yb, II] = sort(yc,'ascend'); % sort the array for processing
            xb = xc(II);
            
            % begin the process to pick other two points from the original 3 point
            % xb(1:3), yb(1:3)
            
            if (y(1) > SOCEndTrg) % y(1)> SOCEndTrg, we need to pick other points < SOCEndTrg
                if ((yb(2) < SOCEndTrg) && (SOCEndTrg < yb(3)))  % yb(2) < SOCEndTrg, pick yb(2)
                    x(2) = xb(2); y(2) = yb(2);
                    if (abs(yb(1)-SOCEndTrg) < abs(yb(3)-SOCEndTrg)) % pick last point from xb(1), and xb(3), according to shortest distance to SOCEndTrg
                        x(3) = xb(1); y(3) = yb(1);
                    else
                        x(3) = xb(3); y(3) = yb(3);
                    end
                end
                if ((yb(1) < SOCEndTrg) && (SOCEndTrg < yb(2))) % yb(1) < SOCEndTrg, pick yb(1)
                    x(2) = xb(1); y(2) = yb(1);
                    if (abs(yb(2)-SOCEndTrg) < abs(yb(3)-SOCEndTrg)) % pick last point from xb(2) and xb(3), according to shortest distance to SOCEndTrg
                        x(3) = xb(2); y(3) = yb(2);
                    else
                        x(3) = xb(3); y(3) = yb(3);
                    end
                end
            end
            if (y(1) < SOCEndTrg) % y(1)< SOCEndTrg, we need to pick other points > SOCEndTrg
                if ((yb(2) < SOCEndTrg) && (SOCEndTrg < yb(3))) % yb(3) > SOCEndTrg, pick yb(3)
                    x(2) = xb(3); y(2) = yb(3);
                    if (abs(yb(1)-SOCEndTrg) < abs(yb(2)-SOCEndTrg)) % pick last point from xb(1) and xb(2), according to shortest distance to SOCEndTrg
                        x(3) = xb(1); y(3) = yb(1);
                    else
                        x(3) = xb(2); y(3) = yb(2);
                    end
                end
                if ((yb(1) < SOCEndTrg) && (SOCEndTrg < yb(2))) % yb(2) > SOCEndTrg, pick this point
                    x(2) = xb(2); y(2) = yb(2);
                    if (abs(yb(1)-SOCEndTrg) < abs(yb(3)-SOCEndTrg)) % pick last point from xb(1) and xb(3), according to shortest distance to SOCEndTrg
                        x(3) = xb(1); y(3) = yb(1);
                    else
                        x(3) = xb(3); y(3) = yb(3);
                    end
                end
            end
        end
        
        y_abs = abs(y-SOCEndTrg);
        [yb, II] = sort(y_abs,'ascend');
        xb = x(II);
        ECMS_s = xb (1);
        
    end
    
    hold off;
    
    toc;
    
    if solution_found
        fprintf ('Search converged. ECMS_s parameter updated in model. \n');
        fprintf ('ECMS_s = %f, SOCEndTrg = %f, SOC_end = %f \n',ECMS_s, SOCEndTrg, SOC_end);
    else
        fprintf ('Search failed to converge. An approximate ECMS_s is updated in the model. \n');
        fprintf ('Refer to the Troubleshooting section of the example page for recommendations. \n');
    end
    
    ECMS_s_tune = ECMS_s;
    % reset model
    Simulink.sdi.markSignalForStreaming(x1,'off');
    load_system(ctrlMdl)
    set_param(ctrlMdl+"/ECMS",'ECMS_s',baseParamName)
    save_system(mdlName,[],'SaveDirtyReferencedModels','on');
    
    % update model and sim
    hws = get_param(ctrlMdl, 'modelworkspace');% get the workspace
    hws.assignin('ECMS_s',ECMS_s);
    save_system(ctrlMdl,[],'SaveDirtyReferencedModels','on');
    open_system(mdlName);
    load_system(battMdl);
    ddObj1 = Simulink.data.dictionary.open('VirtualVehicleTemplate.sldd\BatteryDCDC.sldd');
    dataSectionObj1 = getSection(ddObj1, 'Design Data');
    battChrgMaxValue_obj = getEntry(dataSectionObj1, 'PlntBattChrgMax');
    if  isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj; end
    if ~isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj.getValue; end
    in = in.setVariable('BattCapInit', battChrgMaxValue*SOCinit,'Workspace',battMdl);
    in = in.setVariable('SOCTrgt', SOCEndTrg,'Workspace',ctrlMdl);
    in = in.setVariable('SOCmin', max(SOCEndTrg-20,20.5),'Workspace',ctrlMdl);
    in = in.setVariable('SOCmax', min(SOCEndTrg+20,100),'Workspace',ctrlMdl);
    set_param(ctrlMdl+"/ECMS",'ECMS_s',"ECMS_s");
    save_system(battMdl);
    save_system(ctrlMdl);

end

%%
function in = ModelSetVariable (mdlName, battChrgMaxValue, SOCinit, battMdl, SOCEndTrg, ctrlMdl) 
in = Simulink.SimulationInput(mdlName);
in = in.setVariable('BattCapInit', battChrgMaxValue*SOCinit,'Workspace',battMdl);
in = in.setVariable('SOCTrgt', SOCEndTrg,'Workspace',ctrlMdl);
in = in.setVariable('SOCmin', max(SOCEndTrg-20,20.5),'Workspace',ctrlMdl); 
in = in.setVariable('SOCmax', min(SOCEndTrg+20,100),'Workspace',ctrlMdl);
end

%% 
function [x_out, y_out, ECMS_s, SOC_end, solution_found] = dSOC_generate_starting_points (in, ECMS_CurrentValue, SOCEndTrg)

x_out = zeros (3,1);
y_out = zeros (3,1);
x = zeros (4,1);
y = zeros (4,1);

solution_found = false;
ECMS_s = ECMS_CurrentValue;
alpha = 0.95;

tol = 1;

x1 = ECMS_CurrentValue;  % use current ECMS_s value as the starting point
in = in.setVariable("ECMS_s_tune", x1);
[SOC_end,  dSOC] =dSOCsim_v1(in);  %evaluate x1 results
y1 = SOC_end;
subplot (1,2,1);
plot(x1,dSOC,'bs','MarkerSize',15)
grid on
hold on

subplot (1,2,2);
plot(x1,y1,'bs','MarkerSize',15)
grid on
hold on

if (abs(y1 - SOCEndTrg) <= tol)  % check if a solution found
    ECMS_s  = x1;
    SOC_end = y1;
    solution_found = true;
end

if ((y1 > SOCEndTrg) && ~solution_found)
    x2 = x1 - ECMS_S_dis_up (x1, y1, SOCEndTrg); % use a decreased ECMS_s value as x2
    in = in.setVariable("ECMS_s_tune", x2);
    [SOC_end, dSOC] =dSOCsim_v1(in);
    y2 = SOC_end;

    subplot (1,2,1);
    plot(x2,dSOC,'bs','MarkerSize',15)
    grid on
    hold on

    subplot (1,2,2);
    plot(x2,y2,'bs','MarkerSize',15)
    grid on
    hold on


    if (abs(y2 - SOCEndTrg) <= tol)  % check if a solution found
        ECMS_s  = x2;
        SOC_end = y2;
        solution_found = true;
    end
    if ((y2 > SOCEndTrg) && ~solution_found) % y2 is still too high
        y3_try = SOCEndTrg - 2; % set a lower SOC end trial to it will be easier to get y2 < SOC_EndTrg
        x3 = x1 + alpha*(y3_try - y1)*(x2-x1)/(y2-y1); % use x1, y1, x2, y2 to fit a linear function

$x = x_{1} + α \frac{y - y_{1}}{y_{2} - y_{1}} \times (x_{2} - x_{1})$ such that $x = x_{1} at y = y_{1}, and x = α x_{2} + (1 - α) x_{1} at y = y_{2}$

        dx = ECMS_S_dis_up (x2, y2, SOCEndTrg); % get a pre-defined decrease value
        x3 = min (x3, x2 -   dx);  % upper limit for x3
        x3 = max (x3, x2 - 2*dx);  % lower limit for x3
        in = in.setVariable("ECMS_s_tune", x3); % set the x3 value as ECMS_s
        [SOC_end,  dSOC] =dSOCsim_v1(in);       % evaluate
        y3 = SOC_end;
        subplot (1,2,1);
        plot(x3,dSOC,'bs','MarkerSize',15)
        grid on
        hold on

        subplot (1,2,2);
        plot(x3,y3,'bs','MarkerSize',15)
        grid on
        hold on

        if (abs(y3 - SOCEndTrg) <= tol) % check if a solution is found
        ECMS_s  = x3;
        SOC_end = y3;
        solution_found = true;
        end
        if ~solution_found
            x_out(1) = x1; x_out(2) = x2; x_out(3) = x3; y_out(1) = y1; y_out(2) = y2; y_out(3) = y3; % output 3 points
        end
    end
    if ~solution_found
        if ((y2 > SOCEndTrg) && (y3 >= SOCEndTrg)) % if y3 is still greater than SOCEndTrg, lower it again
            for it_again = 1 : 10
                [x(4)] = Second_order_roots (x_out, y_out, SOCEndTrg - 2);  % decrease again
                dx = ECMS_S_dis_up (x_out (3), y_out (3), SOCEndTrg);   % standard decrease
                x(4) = max (x(4), x_out (3) -   dx);             % make it higher than the standard increase
                x(4) = min (x(4), x_out (3) - 2*dx);             % limit the increase to 2 times the standard
                in = in.setVariable("ECMS_s_tune", x(4));
                [SOC_end, ~] = dSOCsim_v1(in);
                ECMS_s = x(4);
                y(4) = SOC_end;
                x_out (1) = x_out (2);   y_out (1) = y_out (2);
                x_out (2) = x_out (3);   y_out (2) = y_out (3);
                x_out (3) = x(4);        y_out (3) = SOC_end;
                if (abs(y(4) - SOCEndTrg) <= tol); break; end
                if (SOC_end <= SOCEndTrg); break; end
            end
            if (abs(y(4) - SOCEndTrg) <= tol)
                ECMS_s  = x(4);
                SOC_end = y(4);
                solution_found = true;
            end
        end
    end

    if ((y2 < SOCEndTrg) && ~solution_found) % y2 is lower than SOCEndTrg while y1 is greater than SOCEndTrg
        x_min = min (x1, x2); x_max = max (x1, x2); y_min = min(y1, y2); y_max = max (y1, y2); % sort the min and max
        w_min = abs (y_min - SOCEndTrg) / (abs(y_min - SOCEndTrg) + abs(y_max - SOCEndTrg)); % weighting factor
        x3 = (1-w_min) * x_min + w_min * x_max; % x3 is a weighted interpolation between x1 and x2
        in = in.setVariable("ECMS_s_tune", x3);
        [SOC_end,  dSOC] =dSOCsim_v1(in);
        y3 = SOC_end;
        
        subplot (1,2,1);
        plot(x3,dSOC,'bs','MarkerSize',15)
        grid on
        hold on

        subplot (1,2,2);
        plot(x3,y3,'bs','MarkerSize',15)
        grid on
        hold on

        if (abs(y3 - SOCEndTrg) <= tol)
        ECMS_s  = x3;
        SOC_end = y3;
        solution_found = true;
        end
    end
    if ~solution_found
        x_out(1) = x1; x_out(2) = x2; x_out(3) = x3; y_out(1) = y1; y_out(2) = y2; y_out(3) = y3; % output 3 points
    end
end

if ((y1 < SOCEndTrg) && ~solution_found) % y1 < SOCEndTrg
    x2 = x1 + ECMS_S_dis_up (x1, y1, SOCEndTrg); % use a higher ECMS_s value as x2
    in = in.setVariable("ECMS_s_tune", x2);
    [SOC_end, dSOC] =dSOCsim_v1(in);
    y2 = SOC_end;
    subplot (1,2,1);
    plot(x2,dSOC,'bs','MarkerSize',15)
    grid on
    
    subplot (1,2,2);
    plot(x2,y2,'bs','MarkerSize',15)
    grid on
    hold on

    if (abs(y2 - SOCEndTrg) <= tol)
        ECMS_s  = x2;
        SOC_end = y2;
        solution_found = true;
    end
    if ((y2 < SOCEndTrg) && ~solution_found) % x2 is not high enough
        y3_try = SOCEndTrg + 2;                    % set a higher SOC target

$x = x_{1} + α \frac{y - y_{1}}{y_{2} - y_{1}} (x_{2} - x_{1}), x = x_{1}, at y = y_{1}, x = α x_{2} + (1 - α) x_{1} at y = y_{2}$ ,

        x3 = x1 + alpha*(y3_try - y1)*(x2-x1)/(y2-y1); % use x1, y1, x2, y2 as a linear prdeiction
        dx = ECMS_S_dis_up (x2, y2, SOCEndTrg);    % standard increase
        x3 = max (x3, x2 +   dx);                  % lower limit for x3
        x3 = min (x3, x2 + 2*dx);                  % upper limit for x3
        in = in.setVariable("ECMS_s_tune", x3);
        [SOC_end, dSOC] =dSOCsim_v1(in);
        y3 = SOC_end;
        subplot (1,2,1);
        plot(x3,dSOC,'bs','MarkerSize',15)
        grid on

        subplot (1,2,2);
        plot(x3,y3,'bs','MarkerSize',15)
        grid on
        
        if (abs(y3 - SOCEndTrg) <= tol)
            ECMS_s  = x3;
            SOC_end = y3;
            solution_found = true;
        end
        
        if ~solution_found
            x_out(1) = x1; x_out(2) = x2; x_out(3) = x3; y_out(1) = y1; y_out(2) = y2; y_out(3) = y3; % output 3 points
        end
    end
    if ~solution_found
        if ((y2 < SOCEndTrg) && (y3 <= SOCEndTrg))  % y3 is still not high enough
            for it_again = 1 : 10
                [x(4)] = Second_order_roots (x_out, y_out, SOCEndTrg + 2);  % increase again
                dx = ECMS_S_dis_up (x_out (3), y_out (3), SOCEndTrg);   % standard increase
                x(4) = max (x(4), x_out (3) +   dx);             % make it higher than the standard increase
                x(4) = min (x(4), x_out (3) + 2*dx);             % limit the increase to 2 times the standard
                in = in.setVariable("ECMS_s_tune", x(4));
                [SOC_end, dSOC] = dSOCsim_v1(in);
                ECMS_s = x(4);
                y(4) = SOC_end;
                x_out (1) = x_out (2);   y_out (1) = y_out (2);
                x_out (2) = x_out (3);   y_out (2) = y_out (3);
                x_out (3) = x(4);        y_out (3) = SOC_end;
                if (SOC_end >= SOCEndTrg); break; end
                if (abs(y(4) - SOCEndTrg) <= tol); break; end
            end

            subplot (1,2,1);
            plot(x(4),dSOC,'bs','MarkerSize',15)
            grid on

            subplot (1,2,2);
            plot(x(4),y(4),'bs','MarkerSize',15)
            grid on

            if (abs(y(4) - SOCEndTrg) <= tol)
                ECMS_s  = x(4);
                SOC_end = y(4);
                solution_found = true;
            end
        end
    end
    if ((y2 > SOCEndTrg) && ~solution_found)  % y2 is good
        x_min = min (x1, x2); x_max = max (x1, x2); y_min = min(y1, y2); y_max = max (y1, y2); % get min and max for x1, x2, y1, y2
        w_min = abs (y_min - SOCEndTrg) / (abs(y_min - SOCEndTrg) + abs(y_max - SOCEndTrg));   % weighted factor
        x3= (1-w_min) * x_min + w_min * x_max;  % get weighted interpolation between x1, and x2
        in = in.setVariable("ECMS_s_tune", x3);
        [SOC_end, dSOC] =dSOCsim_v1(in);
        y3 = SOC_end;
        subplot (1,2,1);
        plot(x3,dSOC,'bs','MarkerSize',15)
        grid on

        subplot (1,2,2);
        plot(x3,y3,'bs','MarkerSize',15)
        grid on

        if (abs(y3 - SOCEndTrg) <= tol)
            ECMS_s  = x3;
            SOC_end = y3;
            solution_found = true;
        end
    end
    if ~solution_found
        x_out(1) = x1; x_out(2) = x2; x_out(3) = x3; y_out(1) = y1; y_out(2) = y2; y_out(3) = y3; % output three points
    end
end

end
%%
function [SOC_end, dSOC]=dSOCsim_v1(in)

    evalc('out=sim(in)'); %Simulate suppressing build info

    if isempty(out.ErrorMessage)
        SOC=out.logsout.getElement('Battery SOC').Values.Data;
        dSOC=SOC(end)-SOC(1);
        SOC_end = SOC(end);
    else
        dSOC=NaN;
    end
    
end
%%
function dx = ECMS_S_dis_up (ECMS_s, y, SOC_target)

% function computes standard dx = increase/decrease in ECMS_s
% the dx is proportional to the distance between SOC_end (y) and SOC_Target
% also, dx is proportional to ECMS_s, a, y, ff can be adjusted

ff = 0.05;
a  = 0.015;
c = 3.5;
b = ECMS_s / c;   % ECMS_s effect
if (ECMS_s < c-0.1); b = 1; end
dx = a + b * abs(y-SOC_target) * ff;

end
%%
function [z] = Second_order_roots (x, y, y_tar)

$solve the equation y_{t a r} = a z^{2} + b z + c, with y_{i} = a x_{i}^{2} + b x_{i} + c, for 1 \leq i \leq 3 .$

d1 = y(1) / (x(1) - x(2)) / (x(1) - x(3));
d2 = y(2) / (x(2) - x(1)) / (x(2) - x(3));
d3 = y(3) / (x(3) - x(1)) / (x(3) - x(2));

AA = (d1 + d2 + d3);
BB = - (d1 * (x(2) + x(3)) + d2 * (x(1) + x(3)) + d3 * (x(1) + x(2)));
CC = d1 * x(2) * x(3) + d2 * x(1) * x(3) + d3 * x(1) * x(2) - y_tar;

z1 = (-BB + sqrt (BB*BB - 4 * AA * CC)) / 2 / AA;
z = real(z1);

end

Calibrate ECMS Block for Objective Hybrid Vehicle Fuel Economy Assessment

ECMS Calibration Script

Choose HEV Reference Application

ECMS_s Calibration Function

See Also

Related Topics