Calibrate ECMS Block

This example shows how to calibrate the ECMS block in the HEV P2 reference application.

HEV P2 Reference Application

Open the P2 reference application.

autoblkHevP2Start

Set Path

Sets path to the calibration script.

path = fullfile(matlabroot,'examples','autoblks','main');
addpath(path);

Execute Calibration Script

% calibration_ecms_script
pxName="P2"; % PT configuration type, P0,P1,P2,P3,P4
battMdl = "BattHev"+pxName; % battery model - for SOC sweep
mdlName="Hev"+pxName+"ReferenceApplication"; % Hev model
ctrlMdl="Hev"+pxName+"OptimalController"; % Controller model
maxIterat=4;% maximum iterations
SOCinit = 0.85; % initial SOC, unit is in [0, 1]
SOCEndTrg = SOCinit*100;

% 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 + "/" + "Drive Cycle Source";
m = get_param(blk, 'MaskObject');
ECMS_CurrentCycle = m.Parameters(1,1).Value; % Current cycle setting.

mdlWks = get_param(ctrlMdl,'ModelWorkspace');  % find model workspace
% get current ECMS_s value
ECMS_obj = getVariable(mdlWks,'ECMS_s');
if  isnumeric(ECMS_obj); ECMS_CurrentValue = ECMS_obj; end
if ~isnumeric(ECMS_obj); ECMS_CurrentValue = ECMS_obj.Value; end

% extract initial value/name of ECMS tuning parameter
p = Simulink.Mask.get(ctrlMdl+"/ECMS");
baseParamName=p.getParameter("ECMS_s").Value;
battChrgMaxValue_obj = getVariable(get_param(battMdl,'modelworkspace'),'BattChargeMax');
if  isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj; end
if ~isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj.Value; 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')
sgtitle(ECMS_CurrentCycle);
% 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);

### Starting serial model reference simulation build.
### Successfully updated the model reference simulation target for: BattHevP2
### Model reference simulation target for DrivetrainHevP2 is up to date.
### Successfully updated the model reference simulation target for: HevP2OptimalController
### Model reference simulation target for HevP2TransmissionController is up to date.
### Model reference simulation target for MotMappedP2 is up to date.
### Model reference simulation target for SiEngineController is up to date.
### Model reference simulation target for SiMappedEngine is up to date.
### Model reference simulation target for StarterSystemP2 is up to date.

Build Summary

Simulation targets built:

Model                   Action                        Rebuild Reason                                             
=================================================================================================================
BattHevP2               Code generated and compiled.  Change detected in model workspace for model 'BattHevP2'.  
HevP2OptimalController  Code generated and compiled.  Model or library HevP2OptimalController changed.           

2 of 8 models built (6 models already up to date)
Build duration: 0h 3m 2.539s
x1 = 3.430000, y1 = 76.501288 
### Starting serial model reference simulation build.
### Model reference simulation target for BattHevP2 is up to date.
### Model reference simulation target for DrivetrainHevP2 is up to date.
### Model reference simulation target for HevP2OptimalController is up to date.
### Model reference simulation target for HevP2TransmissionController is up to date.
### Model reference simulation target for MotMappedP2 is up to date.
### Model reference simulation target for SiEngineController is up to date.
### Model reference simulation target for SiMappedEngine is up to date.
### Model reference simulation target for StarterSystemP2 is up to date.

Build Summary

0 of 8 models built (8 models already up to date)
Build duration: 0h 0m 26.668s
x2 = 3.861437, y2 = 80.150757 
### Starting serial model reference simulation build.
### Model reference simulation target for BattHevP2 is up to date.
### Model reference simulation target for DrivetrainHevP2 is up to date.
### Model reference simulation target for HevP2OptimalController is up to date.
### Model reference simulation target for HevP2TransmissionController is up to date.
### Model reference simulation target for MotMappedP2 is up to date.
### Model reference simulation target for SiEngineController is up to date.
### Model reference simulation target for SiMappedEngine is up to date.
### Model reference simulation target for StarterSystemP2 is up to date.

Build Summary

0 of 8 models built (8 models already up to date)
Build duration: 0h 0m 25.533s
x3 = 4.422918, y3 = 82.791282 
### Starting serial model reference simulation build.
### Model reference simulation target for BattHevP2 is up to date.
### Model reference simulation target for DrivetrainHevP2 is up to date.
### Model reference simulation target for HevP2OptimalController is up to date.
### Model reference simulation target for HevP2TransmissionController is up to date.
### Model reference simulation target for MotMappedP2 is up to date.
### Model reference simulation target for SiEngineController is up to date.
### Model reference simulation target for SiMappedEngine is up to date.
### Model reference simulation target for StarterSystemP2 is up to date.

Build Summary

0 of 8 models built (8 models already up to date)
Build duration: 0h 0m 23.503s
x(4) = 4.732032, y(4) = 84.960746 

if (solution_found == 0)
    
    for i = 1 : maxIterat
        
        [yy, II] = sort(y,'ascend'); % sort the x, y array for plotting
        xx = x(II);
        subplot (1,2,2)
        if (i == 1); plot(xx,yy,'b--o','LineWidth',2); end
        if (i == 2); plot(xx,yy,'g--o','LineWidth',2); end
        if (i == 3); plot(xx,yy,'r--o','LineWidth',2); end
        if (i == 4); plot(xx,yy,'y--o','LineWidth',2); end
        xlabel('ECMS\_s')
        ylabel('SOC')
        hold off
        % solve second order equation to get z = ECMS_s corresponds to y = SOCEndTrg
        

$$\begin{array}{l}SOCEndTrg=SO{C}_{init}=a{x}^2+bx+c,\text{where}a,b,c\text{satisfies}{y}_i=a{x}_i^2+b{x}_i+c,1\le i\le 3.\text{with}\\ x_1={\text{ECMS}}_s^{(1)},y_1=SO{C}_{end}({\text{ECMS}}_s^{(1)}),x_2={\text{ECMS}}_s^{(2)},y_2=SO{C}_{end}({\text{ECMS}}_s^{(2)}),x_3={\text{ECMS}}_s^{(3)},y_3=SO{C}_{end}({\text{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)
        hold on
        fprintf ('i= %u, ECMS_s = %f, SOC_end = %f \n',i, z, SOC_end);
        fprintf ('i= %u, x1 = %f, x2 = %f, x3 = %f, \n',i, x(1), x(2), x(3));
        fprintf ('i= %u, y1 = %f, y2 = %f, y3 = %f, \n',i, y(1), y(2), y(3));
        
        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;

Elapsed time is 656.183812 seconds.

if (solution_found == 1)
    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);
end

Search converged. ECMS_s parameter updated in model. 

ECMS_s = 4.732032, SOCEndTrg = 85.000000, SOC_end = 84.960746 

if (solution_found == 0)
    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);
battChrgMaxValue_obj = getVariable(get_param(battMdl,'modelworkspace'),'BattChargeMax');
if  isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj; end
if ~isnumeric(battChrgMaxValue_obj); battChrgMaxValue = battChrgMaxValue_obj.Value; 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);
% sim(in);
% open_system(mdlName+"/Visualization/Performance and FE Scope");

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 = 0;
ECMS_s = ECMS_CurrentValue;
alpha = 0.8;

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)
hold on
fprintf ('x1 = %f, y1 = %f \n',x1, y1);
if (abs(y1 - SOCEndTrg) <= tol)  % check if a solution found
    ECMS_s  = x1;
    SOC_end = y1;
    solution_found = 1;
end

if ((y1 > SOCEndTrg) && (solution_found == 0))
    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)
    hold on
    fprintf ('x2 = %f, y2 = %f \n',x2, y2);
    if (abs(y2 - SOCEndTrg) <= tol)  % check if a solution found
        ECMS_s  = x2;
        SOC_end = y2;
        solution_found = 1;
    end
    if ((y2 > SOCEndTrg) && (solution_found == 0)) % 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+\alpha \frac{y-y_1}{y_2-y_1}\times (x_2-x_1)$$ such that $$x=x_1\text{at}y=y_1,\text{and}x=\alpha x_2+(1-\alpha )x_1\text{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)
        hold on
        fprintf ('x3 = %f, y3 = %f \n',x3, y3);
        if (abs(y3 - SOCEndTrg) <= tol) % check if a solution is found
        ECMS_s  = x3;
        SOC_end = y3;
        solution_found = 1;
        end
        if (solution_found == 0)
            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 == 0)
        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 = 1;
            end
        end
    end

    if ((y2 < SOCEndTrg) && (solution_found == 0)) % 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)
        hold on
        fprintf ('x3 = %f, y3 = %f \n',x3, y3);
        if (abs(y3 - SOCEndTrg) <= tol)
        ECMS_s  = x3;
        SOC_end = y3;
        solution_found = 1;
        end
    end
    if (solution_found == 0)
        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 == 0)) % 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)
    fprintf ('x2 = %f, y2 = %f \n',x2, y2);
    if (abs(y2 - SOCEndTrg) <= tol)
        ECMS_s  = x2;
        SOC_end = y2;
        solution_found = 1;
    end
    if ((y2 < SOCEndTrg) && (solution_found == 0)) % x2 is not high enough
        y3_try = SOCEndTrg + 2;                    % set a higher SOC target
         

$$x=x_1+\alpha \frac{y-y_1}{y_2-y_1}(x_2-x_1),x=x_1,\text{at}y=y_1,x=\alpha x_2+(1-\alpha )x_1\text{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)
        fprintf ('x3 = %f, y3 = %f \n',x3, y3);
        if (abs(y3 - SOCEndTrg) <= tol)
            ECMS_s  = x3;
            SOC_end = y3;
            solution_found = 1;
        end
        if (solution_found == 0)
            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 == 0)
        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)
            fprintf ('x(4) = %f, y(4) = %f \n',x(4), y(4));
            if (abs(y(4) - SOCEndTrg) <= tol)
                ECMS_s  = x(4);
                SOC_end = y(4);
                solution_found = 1;
            end
        end
    end
    if ((y2 > SOCEndTrg) && (solution_found == 0))  % 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)
        fprintf ('x3 = %f, y3 = %f \n',x3, y3);
        if (abs(y3 - SOCEndTrg) <= tol)
            ECMS_s  = x3;
            SOC_end = y3;
            solution_found = 1;
        end
    end
    if (solution_found == 0)
        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)

out=sim(in);
% for ii=1:length(out)
    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, c, 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)

$$\text{solve}\text{the}\text{equation}{y}_{tar}=a{z}^2+bz+c,\text{with}{y}_i=a{x}_i^2+b{x}_i+c,\text{for}1\le i\le 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