Problem with constraints SoC - battery storage system
Mislav Hruškar
Dear John,
I am making an MPC of battery energy storage for application in railway and my dynamics are PWA. I have 10 region - 10 PWA models wich describe normalized
state-of-charge (SoC) – from the interval [ 0.1 – 1] .
PWA modell goes to region 1 if -x1_Pc <=u{k}(2,1)<= - 1e-12), to region 2 if -x2_Pc
<=u{k}(2,1)<= -x1_Pc - 1e-12) and so on. u{k} (1,1) is discharging power component of storage system and u{k} (2,1) is charging power component. Charging and discharging cannot be at the same time and only one model of PWA can be active at any time. d{k} is electricity price for time instant k, and r{k} denotes power consumption or production.
Now, the problem is why although I bound variable x{k} and when I spin MPC algortihm through the 300 steps in 5th step x{5} are 0.05 and u{5} = NaN. Also u{6}….u{300} are NaN and x{5}…x{300} are 0.05.
I made print of x in 1st step of MPC and I saw that x which I sent to optimization (initial state) is not the same (its zero) and other x are 0.1. u{k} (1,1) are max value although it should be 0 because the battery is discharged after a few steps.
I do not know what to do anymore, I lost a lot of time on this, and therefore I am asking for help.
Thank you in advance for any help you can provide.
Best regards,
Mislav
Here is the code:
% Model data
A = 1;
B1 = (delta_T/E0)*[0, -theta_Pc(index_Pc_(1),1), 0];
B2 = (delta_T/E0)*[0, -theta_Pc(index_Pc_(2),1), 0];
B3 = (delta_T/E0)*[0, -theta_Pc(index_Pc_(3),1), 0];
B4 = (delta_T/E0)*[0, -theta_Pc(index_Pc_(4),1), 0];
B5 = (delta_T/E0)*[0, -theta_Pc(index_Pc_(5),1), 0];
B6 = (delta_T/E0)*[-theta_Pd(index_Pd_(1),1), 0, 0];
B7 = (delta_T/E0)*[-theta_Pd(index_Pd_(2),1), 0, 0];
B8 = (delta_T/E0)*[-theta_Pd(index_Pd_(3),1), 0, 0];
B9 = (delta_T/E0)*[-theta_Pd(index_Pd_(4),1), 0, 0];
B10 = (delta_T/E0)*[-theta_Pd(index_Pd_(5),1), 0, 0];
b1 = (delta_T/E0)*[theta_Pc(index_Pc_(1),2)];
b2 = (delta_T/E0)*[theta_Pc(index_Pc_(2),2)];
b3 = (delta_T/E0)*[theta_Pc(index_Pc_(3),2)];
b4 = (delta_T/E0)*[theta_Pc(index_Pc_(4),2)];
b5 = (delta_T/E0)*[theta_Pc(index_Pc_(5),2)];
b6 = (delta_T/E0)*[theta_Pd(index_Pd_(1),2)];
b7 = (delta_T/E0)*[theta_Pd(index_Pd_(2),2)];
b8 = (delta_T/E0)*[theta_Pd(index_Pd_(3),2)];
b9 = (delta_T/E0)*[theta_Pd(index_Pd_(4),2)];
b10 = (delta_T/E0)*[theta_Pd(index_Pd_(5),2)];
nx = 1; % Number of states
nu = 3; % Number of inputs
% MPC data
load energy_price_profil.mat
load Ptr_5h.mat
N = 300;
figure(3)
stairs(price(:,:),'b');
grid on;
hold on;
stairs(PTR(:,:), 'r');
legend('Energy price', 'Ptr');
xlabel('k');
% Initial state
x0 = 0.7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u = sdpvar(repmat(nu,1,N),repmat(1,1,N));
x = sdpvar(repmat(nx,1,N+1),repmat(1,1,N+1));
r = sdpvar(1,N);
d = sdpvar(1,N);
mod = binvar(repmat(10,1,N),ones(1,N));
uu = binvar(repmat(10,1,N),repmat(1,1,N));
ucase = binvar(repmat(2,1,N),repmat(1,1,N));
xcase = binvar(repmat(1,1,N),repmat(1,1,N));
constraints = [];
objective = 0;
for k = 1:N
XX = [' Ovo je postavljanje optimizacijskog problema, korak: ',num2str(k)];
disp(XX);
objective = objective + d(k)*( r(k) + [-1, -1, 1]*u{k} )*delta_T;
constraints = [constraints, implies(mod{k}(1), x{k+1} == A*x{k} + B1*u{k} + b1)];
constraints = [constraints, implies(mod{k}(2), x{k+1} == A*x{k} + B2*u{k} + b2)];
constraints = [constraints, implies(mod{k}(3), x{k+1} == A*x{k} + B3*u{k} + b3)];
constraints = [constraints, implies(mod{k}(4), x{k+1} == A*x{k} + B4*u{k} + b4)];
constraints = [constraints, implies(mod{k}(5), x{k+1} == A*x{k} + B5*u{k} + b5)];
constraints = [constraints, implies(mod{k}(6), x{k+1} == A*x{k} + B6*u{k} + b6)];
constraints = [constraints, implies(mod{k}(7), x{k+1} == A*x{k} + B7*u{k} + b7)];
constraints = [constraints, implies(mod{k}(8), x{k+1} == A*x{k} + B8*u{k} + b8)];
constraints = [constraints, implies(mod{k}(9), x{k+1} == A*x{k} + B9*u{k} + b9)];
constraints = [constraints, implies(mod{k}(10), x{k+1} == A*x{k} + B10*u{k} + b10)];
constraints = [constraints, implies(uu{k}(1), -x1_Pc <=u{k}(2,1)<= - 1e-12)];
constraints = [constraints, implies(uu{k}(2), -x2_Pc <=u{k}(2,1)<= -x1_Pc - 1e-12) ];
constraints = [constraints, implies(uu{k}(3), -x3_Pc <=u{k}(2,1)<= -x2_Pc - 1e-12)];
constraints = [constraints, implies(uu{k}(4), -x4_Pc <=u{k}(2,1)<= -x3_Pc - 1e-12)];
constraints = [constraints, implies(uu{k}(5), -Pc_max <=u{k}(2,1)<= -x4_Pc - 1e-12)];
constraints = [constraints, implies(uu{k}(6), 0 <=u{k}(1,1)<= x1_Pd - 1e-12)];
constraints = [constraints, implies(uu{k}(7), x1_Pd <=u{k}(1,1)<= x2_Pd - 1e-12)];
constraints = [constraints, implies(uu{k}(8), x2_Pd <=u{k}(1,1)<= x3_Pd - 1e-12)];
constraints = [constraints, implies(uu{k}(9), x3_Pd <=u{k}(1,1)<= x4_Pd - 1e-12)];
constraints = [constraints, implies(uu{k}(10), x4_Pd <=u{k}(1,1)<= Pd_max)];
constraints = [constraints, implies(ucase{k}(1), u{k}(1,1)==0)];
constraints = [constraints, implies(ucase{k}(2), u{k}(2,1)==0)];
constraints = [constraints, implies(xcase{k}(1), 0.1<=x{k}(:)<=1)];
constraints = [constraints, mod{k}(1) + mod{k}(2) + mod{k}(3) + mod{k}(4) + mod{k}(5) == ucase{k}(1)];
constraints = [constraints, mod{k}(6) + mod{k}(7) + mod{k}(8) + mod{k}(9) + mod{k}(10) == ucase{k}(2)];
constraints = [constraints, xcase{k}(1) == 1];
constraints = [constraints, mod{k}(1) == uu{k}(1)];
constraints = [constraints, mod{k}(2) == uu{k}(2)];
constraints = [constraints, mod{k}(3) == uu{k}(3)];
constraints = [constraints, mod{k}(4) == uu{k}(4)];
constraints = [constraints, mod{k}(5) == uu{k}(5)];
constraints = [constraints, mod{k}(6) == uu{k}(6)];
constraints = [constraints, mod{k}(7) == uu{k}(7)];
constraints = [constraints, mod{k}(8) == uu{k}(8)];
constraints = [constraints, mod{k}(9) == uu{k}(9)];
constraints = [constraints, mod{k}(10) == uu{k}(10)];
constraints = [constraints, [0; -Pc_max; 0]<=u{k}<=[Pc_max; 0; 16.5], -3<=(r(k) + [-1, -1, 1]*u{k})<=16.5, 0.1<=x{k}<=1];
end
constraints = [constraints, 0.1<=x{N+1}<=1];
parameters_in = {x{1},r,d};
solutions_out = {[u{:}], [x{:}]};
options = sdpsettings('verbose',1,'solver','cplex','cplex.qpmethod',1,'cachesolvers',1);
controller = optimizer(constraints, objective, options, parameters_in, solutions_out);
x = x0;
U_=[];
x_=x0;
for i = 1:1
rr = PTR(i:i+N-1);
dd = price(i:i+N-1);
[solutions, diagnostics] = controller{{x, rr, dd}};
U = solutions{1};
X = solutions{2};
U_=[U_,U(:,1)];
if (-x1_Pc<=U(2,1)&&U(2,1)<0)
pr1 = [' Korak: ', num2str(i), 'Model 1 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr1);
x = A*x + B1*U(:,1) + b1;
elseif (-x2_Pc<=U(2,1)&&U(2,1)<-x1_Pc)
pr2 = [' Korak: ', num2str(i), 'Model 2 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr2);
x = A*x + B2*U(:,1) + b2;
elseif (-x3_Pc <=U(2,1)&&U(2,1)< -x2_Pc)
pr3 = [' Korak: ', num2str(i), 'Model 3 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr3);
x = A*x + B3*U(:,1) + b3;
elseif (-x4_Pc <=U(2,1)&&U(2,1)< -x3_Pc)
pr4 = [' Korak: ', num2str(i), 'Model 4 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr4);
x = A*x + B4*U(:,1) + b4;
elseif (-Pc_max <=U(2,1)&&U(2,1)< -x4_Pc)
pr5 = [' Korak: ', num2str(i), 'Model 5 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr5);
x = A*x + B5*U(:,1) + b5;
elseif (0 <U(1,1)&&U(1,1)< x1_Pd)
pr6 = [' Korak: ', num2str(i), 'Model 6 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr6);
x = A*x + B6*U(:,1) + b6;
elseif (x1_Pd <=U(1,1)&&U(1,1)< x2_Pd)
pr7 = [' Korak: ', num2str(i), 'Model 7 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr7);
x = A*x + B7*U(:,1) + b7;
elseif (x2_Pd <=U(1,1)&&U(1,1)< x3_Pd)
pr8 = [' Korak: ', num2str(i), 'Model 8 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr8);
x = A*x + B8*U(:,1) + b8;
elseif (x3_Pd <=U(1,1)&&U(1,1)< x4_Pd)
pr9 = [' Korak: ', num2str(i), 'Model 9 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr9);
x = A*x + B9*U(:,1) + b9;
elseif (x4_Pd <=U(1,1)&&U(1,1)<=Pd_max)
pr10 = [' Korak: ', num2str(i), 'Model 10 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr10);
x = A*x + B10*U(:,1) + b10;
elseif (0==U(1,1)&&U(2,1)==0)
pr11 = [' Korak: ', num2str(i), 'Model 11 U(1,1)=', num2str(U(1,1)), 'U(2,1)=', num2str(U(2,1))];
disp(pr11);
x = A*x ;
end
x_=[x_,x];
end
Johan Löfberg
Mislav Hruškar
Hi Johan,
thank you very much for your speedy reply.
I have put my MATLAB code in a .ZIP file , and all you need is to run the m-file MPC_Yalmip. I’ve set the prediction horizon to N = 300 steps. The number of steps for which solves the problem of minimum of cost function is also 300 steps. You will see that after 13th step state variable x has a value below 0.1. (between 1st and 13th we have discharging - decreasing of x), in 14th step x=0.0772, although it should not reach below 0.1 because that value is set on low limit (constraint x is from interval [0.1, 1]). After 13th step value x should be increasing or stable, definitely not falling.
If you change in 163rd line: for i = 1:1 (just one cycle) , run code and observe variable x (I made it so it can be seen - solutions_out = {[u{:}], [x{:}]};) , you we will see that x{1}(1)=0, x{2}(1)...x(300)(1)=0.1. I don't now why, but I have set initial state x{1}(1) = x0{1}(1) = 0.7 which is sent to solver, and in no way it should be 0. I tried everything and don’t know what to do more, so please help me.
Best regards and many thanks,
Mislav Hruskar
Johan Löfberg
Johan Löfberg
Mislav Hruškar
yes, the problem was that I didn't use constraint: constraints = [constraints, sum(mod{k})==1];.
Now I put this one condition and I changed cost function. Also I added in front of expression abs so it all works perfectly.
Thank you very much for your prompt reply and best regards,
Mislav Hruskar
