Robust MPC infeasible from Dr. Lofberg's example for extreme realizations of w = 1 for all t

[Alex Keller]

Hello Professor,

I am attempting to learn your RMPC closed-loop approximation method and have run into a problem which I cannot discern the cause (whether my understanding of the optimization problem or a problem with the method). The problem is when simulating your example from /robustmpc/ with the disturbance realization (w=1), the problem becomes infeasible after 4 steps. 

My current understanding is that this is due to the first max_w of J which results in preparing for the worst case disturbance realization of w = -1 since this maximizes the cost function. Thus, the disturbance w=1 is never optimized for. However, this contradicts the fact that the optimization is infeasible at the 4th step (before crossing Y=1), because w=-1 for all t should remain feasible (ie Y should remain below 1).

Could you please help me understand why your example fails with a disturbance within the provided uncertainty bounds. Thank you.

For comparison, the open-loop RMPC response has been attached as well. Note: I have read your IEEE paper where the example comes from as well as your dissertation, but am still a new student to optimization and MPC.

Code for reproducing problem:

%Example Closed Loop RMPC from https://yalmip.github.io/example/robustmpc/

%% NOTE: Lofberg's optimization fails when we realize the worst case disturbances of w = 1 every step

clear all

addpath(genpath('YALMIP-master'))

yalmip('clear')

A = [2.938 -0.7345 0.25;4 0 0;0 1 0];

B = [0.25;0;0];

C = [-0.2072 0.04141 0.07256];

E = [0.0625;0;0];

N = 10;

W = sdpvar(N,1);

x = sdpvar(3,1);

G = [-1 <= W <= 1];

%% Closed loop approximation

V = sdpvar(N,1);

L = sdpvar(N,N,'full').*(tril(ones(N))-eye(N));

U = L*W + V;

Y = [];

xk = x;

for k = 1:N

 xk = A*xk + B*U(k)+E*W(k);

 Y = [Y;C*xk];

end

F = [Y <= 1, -1 <= U <= 1];

objective = norm(Y-1,1) + norm(U,1)*0.01;

[Frobust,h] = robustify([F, G],objective,[],W);

xk = [0;0;0];

ops = sdpsettings;

for i = 1:25

    optimize([Frobust, x == xk(:,end)],h,ops);

    %Simulate with random disturbance same as Lofberg's example

%     xk = [xk A*xk(:,end) + B*value(U(1)) + E*(-1+2*rand(1))];

    %Simulate with worst case disturbance w=1

    xk = [xk A*xk(:,end) + B*value(U(1)) + E*1];

end

figure(1)

hold on

plot(C*xk,'b')

xk = [0;0;0];

ops = sdpsettings;

for i = 1:25

    optimize([Frobust, x == xk(:,end)],h,ops);

    %Simulate with random disturbance same as Lofberg's example

%     xk = [xk A*xk(:,end) + B*value(U(1)) + E*(-1+2*rand(1))];

    %Simulate with worst case disturbance w=-1

    xk = [xk A*xk(:,end) + B*value(U(1)) + E*(-1)];

end

figure(1)

hold on

plot(C*xk,'r')

plot([0,30],[1,1],'g')

legend('Disturbance realizations: +1 for all t','Disturbance realizations: -1 for all t','Constraint')

[Johan Löfberg]

there is nothing in the method which ensures recursive feasibility, it is entirely focused on solving the given minmax poblem. For stability and recursive feasibility much more is needed (or simple standard stuff such as slacks)

[Alex Keller]

Thank you, that clears it up nicely. 

One follow up question: Is there a way to extract the extreme value of W which maximizes the cost function (for some time t when we run the robust mpc with actual state values)?

[Johan Löfberg]

..and it can also be due to numerical issues. The solution returned at iteration k is epsilon infeasible, and then that causes issues in the next iteration when the new solution possible has to be more than epsilon infeasible. Same issue as in any MPC algorithm, uncertainty or not

No, extracting W is not possible, as there might not even be such a value

-t <=W <=t for all |W|<=1 does not have *a* worst-case, it has two completely different cases, and that generalizes