How to implement optimal control using Pontryagin's Minimum Principle for HESS in EV application as in attached paper?

I am trying to implement a paper based on optimal control of energy management of HESS in EV. It's based on Pontryagin's Minimum Principle using Hamiltonian, state and costate equations. How to simulate this? Is it better to use matlab m file or simulink?

Please find the paper attached.

Answer 1

[Stephan]

Hi,

there are some examples here on answers, where people used Matlab with Symbolic Math Toolbox to solve this kind of problems. I do not remember Simulink modesl in this context - but maybe there are some. A newer example using Matlab with a working code you can find here:

https://de.mathworks.com/matlabcentral/answers/438910-how-to-solve-in-matlab-2018b

Best regards

Stephan

[Praveen Kumar]

Thank you. I will go through it and try to execute.

[Praveen Kumar]

clear;
%Initializations and values of constant variables
C = 28;
Rs = 0.01;
Rp = 10e3;
Vsmax = 305;
Preq = 50e3;
k = 1;
Vb = 360;
Rb = 0.3;
p1 = -3;
tf = 5;
% State equations
syms x1 Pb p1; %Pb control variable & x1 = SOC_s is state variable
Dx1 = -1/C*((1/(2*Rs)+1/Rp)*x1-1/(2*Rs)*sqrt(x1^2-4*Rs/Vsmax^2*(Preq-Pb)));
% Cost function inside integral
syms g;
g = k*((Vb - sqrt(Vb^2 - 4*Rb*Pb))/(2*Rb))^2;
% Hamiltonian
syms p1 H;
H = g + p1*Dx1;
% Costate equations
% Dx1 = diff(H,p1);
Dp1 = -diff(H,x1);
% Solve for control Pb
dPb = diff(H,Pb);
sol_Pb = solve(dPb,Pb);
% Substitute Pb to state equations
Dx1 = subs(Dx1,Pb,sol_Pb);
% Convert symbolic objects to strings for using 'dsolve'
eq1 = strcat('Dx1=',char(Dx1));
eq2 = strcat('Dp1=',char(Dp1));
sol_h = dsolve(eq1,eq2);
% Use boundary conditions to determine the coefficients
% Case a: (a)Dx1(0)=90 Dx1(tf)=50
conA1 = 'x1(0) = 0';
conA2 = 'x1(10) = 1';
sol_a = dsolve(eq1,eq2,conA1,conA2);
% plot solutions
figure(1)
ezplot(sol_a.x1,[0 1]); hold on;
ezplot(sol_a.p1,[0 1]); % plot the control u=-p1
axis([0 20 0 1]);
xlabel('time');
ylabel('states');
title('Solution of PMP')

I have tried to implement my example using the framework given in the above link. I am getting an error for dsolve.

Error using mupadengine/feval (line 187)

Invalid number of arguments.

Error in dsolve>mupadDsolve (line 340)

T = feval(symengine,'symobj::dsolve',sys,x,options);

Error in dsolve (line 194)

sol = mupadDsolve(args, options);

Error in PMP_18012019 (line 33)

sol_h = dsolve(eq1,eq2);.

[abs mhr]

hi

were you successful to resolve the errors?

[NADA RIFAI]

Hello,

I have the same type of problem but with constraints, how can I include them in the solution?

Can you please help me solve this problem.

thank you

[Praveen Kumar]

Hi,

I am not exactly sure about including constraints. You may try to avoid those control variable candidates which violates the constraints i.e use only arg min(H) than dH=0; while maintaining the constraints. In some cases, you might have to use jump conditions.

Please search papers on " jump conditions of PMP " which are useful during state constraints.

Also you may see the research works of Namwook Kim (S.Korea, N.Kim PMP) or Simona Onori (Stanford) for jump conditions.

[NADA RIFAI]

Hi,

My problem is similar to yours, except that I have inequality constraints.

You can find the problem in the attached document.

Thank you very much.

[Praveen Kumar]

Yea, I have seen Hemi's work. He has implemented in Simulink while I used m-code. I think these inequality constrains can be implemented in sameway I discussed in last answer.

Suppose, when you calculate the state variables for next instant or at any instant based on the control variable, if it violates any constraints you may use the next higher or lower value of control variable.

Eg: For a particular instant, from minimization of H you found that the optimal battery current Ib = 8A, but if you discharge with 8A then your SOC will go below 0.2 (lower limit). So in such case instead of 8A choose another value maybe 6 or 7A so that SOC limits are not violated. So choosing that current from the available control candidates must be justified.

Hope this helps, as this solved some of the problems. Even I am still working on the problem with constraints, but not solved it efficiently.

You may contact Hemi via Linkedin as he is very helpful (he sent me his thesis) or drop me an email so that I can send his thesis (but in French).

[NADA RIFAI]

Thank you so much. Yes French it's perfect for me. [e-mail]

[LORENZO MARTOCCIA]

Can you send this document (thesis file) to my mail? [e-mail]

Thank you so much.