time varying input ...

[LP]

Hi Johan,

I am back with one more question. You solved my problem last time regarding the understanding of the closed-loop Robust MPC.

The question now is that when using closed loop MPC, the input U is not linear anymore, and now I would like to know, how to construct objective function  of the form:

objective = sum(U(1:N).*var1(1:N)) - sum(R(1:N).*var2(1:N))

in which, var1 and var2 are time varying. As the expression is polynomial, robustify is not able to cope with it. But it's a pity since I already know the exact value of var1, and var2, and the objective funciton should remain in theory bilinear. But I am not able to find a way to call the model as in to roll var1, and var2 with time, such that:

optimizer(constraints, objective,ops,{blah blah, var1, var2},{U,R});

As a reference, following is the code in which var1 = usep, and var2 = res:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

function [Ureal,Rreal] = closedloopRobustMPC_optimizer( set_back, Q_real, Q_mean_mat, Q_mean_vec, Q_ub, Q_lb,N,sim_Time,m,Ad,Bd,usep,res)
usep_in = sdpvar(N,1);
res_in = sdpvar(N,1);

slack = sdpvar(N,1);
W = sdpvar(N,1);
x1 = sdpvar(5,1);
x2 = sdpvar(5,1);
Q_mean_mat_var = sdpvar(4,N);
sb = sdpvar(N,1);
C = [0 0 0 0 1];
Y1 = [];
Y2 = [];
xk1 = x1;
xk2 = x2;
Q_mean_mat2 = Q_mean_mat_var;
set_back2 = sb;

V = sdpvar(N,1);
Z = sdpvar(N,1);
memory = 48; % Between N and 0
L = sdpvar(N,N,'full').*triu(tril(ones(N),-1),-memory);
M= sdpvar(N,N,'full').*triu(tril(ones(N),-1),-memory);
slack = sdpvar(N,1);



U = L*W + V;
R = M*W + Z;



objective= 0;
for k = 1:N
    xk2 = Ad*xk1 - Bd(:,1)*(U(k)-R(k)) + Bd(:,2:end-1)*Q_mean_mat2(:,k)  + Bd(:,end)*W(k);
    Y2 = [Y2;C*xk2];
   
    xk1 = Ad*xk1 - Bd(:,1)*(U(k)) + Bd(:,2:end-1)*Q_mean_mat2(:,k)  + Bd(:,end)*W(k);
    Y1 = [Y1;C*xk1];


end

%The input is bounded, and the objective is to stay as close as possible to
%the level y=1 while guaranteeing that y<1 is satisfied despite the
%disturbances w.
% F = [70 - slack - set_back2 <= Y <= 74 + slack + set_back2, 0 <= U <= 100, slack >= 0];
F = [70 - slack - set_back2 <= Y1 <= 74 + slack + set_back2,70 - slack - set_back2 <= Y2 <= 74 + slack + set_back2, ...
    0 <= U + R <= 100, slack >= 0, R >= 0, U - R >= 0];
% objective = sum(1e6.*slack) + sum(U.*usep(1:N)')*1800*0.00055;
objective = 1e6*sum(slack) + sum(U.*usep_in')*1800*0.00055 - sum(R.*res_in')*1800*0.00055;

ops = sdpsettings('verbose',2, 'solver', '+cplex');

xk1 = [72.896;75.9814;75.8665;72.4428;72.5114];
xk2 = [72.896;75.9814;75.8665;72.4428;72.5114];
set_back2 = set_back;
Q_mean_mat2 = Q_mean_mat;
Ureal = [];
Rreal = [];

G = [uncertain(W), min(Q_lb(5,1:N)) <= W <= max(Q_ub(5,1:N))];

controller = optimizer([F,G], objective,ops,{x1(:,1),x2(:,1), Q_mean_mat_var,sb,usep_in,res_in},{U,R});

for i = 1:48
   

    [sol,diagnostics] = controller{{xk1(:,end), xk2(:,end), Q_mean_mat2(:,i:i+N-1),set_back2(i:i+N-1)',usep(1:N)',res(1:N)'}};
    uk = sol{1}(1,1);
    rk = sol{2}(1,1);
   
   
    xk2 = [xk2 Ad*xk1(:,end) -  Bd(:,1)*(uk-rk) + Bd(:,2:end-1)*Q_mean_mat2(:,i) + Bd(:,end)*Q_real(i)];
    xk1 = [xk1 Ad*xk1(:,end) -  Bd(:,1)*uk + Bd(:,2:end-1)*Q_mean_mat2(:,i) + Bd(:,end)*Q_real(i)];
   
   
    Ureal = [Ureal;uk];
    Rreal = [Rreal;rk];


end
figure();
plot(C*xk1(:,2:end),'r','LineWidth',3.5);  hold on
plot(C*xk2(:,2:end),'k','LineWidth',3.5);
plot(1:sim_Time,74 + set_back2(1:sim_Time)','LineWidth',4);
plot(1:sim_Time,70 - set_back2(1:sim_Time)','LineWidth',4);
title('CL-RMPC');
figure();
stairs(Ureal,'r','LineWidth',3.5); hold on;
stairs(Rreal,'k','LineWidth',3.5);
title('CL-RMPC')


end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Johan Löfberg]

Are you using the latest version of YALMIP. Cannot run your code, but similar works

sdpvar x w c
 P = optimizer([-2 <= x <= 2, uncertain(w), -1 <= w <= 1], x + x*w*c,sdpsettings('solver','+gurobi'),c,x)
Optimizer object with 1 inputs and 1 outputs. Solver: GUROBI-GUROBI
K>> P{0}

ans =

    -2

K>> P{-10}

ans =

     0

[LP]

yalmip('version') =  20141218

The polynomial expression given by you above, even though the principle is exactly the same for my model. But for me, It just get stuck in robustification module and doesn't get out. It doesn't say any error or anything. I have also lowered the parameter variables non zeros entries etc. still no luck.... I don't have gurobi, I use cplex, but I don't think it could be the problem ....

Following is what I get, and it stays like this if I will even leave it for 10 hours.

***** Starting YALMIP robustification module. *********************
 - Detected 48 uncertain variables
 - Detected 48 independent group(s) of uncertain variables
 - Eliminating uncertainty using explicit maximization of inf-norm

[LP]

Here are the files that should reproduce the problem I am facing.

Let me know if you'd need anything.

I still think it is a polynomial expression over 48 time steps which grows huge in robustification module.

Maybe I try doing this manually, if it'll help ...

Do let me know if you could figure out the issue...

Thanks!

[Johan Löfberg]

Start by installing the latest version of yalmip

[Johan Löfberg]

OK, now I see what is happening. Due to the complicated structure with parametric vbariables, YALMIP resorts to an enumeration scheme, which absolutely won't work in this high-dimensional space. At the moment, there is no way to solve this using optimizer

[LP]

Hi Johan, I am back!

Ok, since unable to use "robustify". I set out to solve this manually using bounds and applying strong duality concept.

But unfortunately, I am not able to formulate the parametrized matrix "L" properly in the feedback U = L*W + V.
And, all I get are zeros in the parameter matrix. I hope you can give a better insight into the problem. It will be highly appreciated.

Thanks a lot,

%%%%%%%%%% SYSTEM %%%%%%%%%%%%%%
so the system with respect to the structure of MPC is the same as the one given in the example
in short, this is the system:
obj = min(U) + min(rho*slack) ...
slack >= 0
Sx <= s_x
Su <= s_u

U = L*W + V

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%% Code %%%%%%%%%%%%%%
V = sdpvar(N,1);

memory = 48; % Between N and 0
L = sdpvar(N,N,'full').*triu(tril(ones(N),-1),-memory);

slack = sdpvar(N,1);



lambda1 = sdpvar(5,N);
lambda2 = sdpvar(5,N);
lambda3 = sdpvar(N,N);
lambda4 = sdpvar(N,N);
lambda5 = sdpvar(N,N);
lambda6 = sdpvar(N,N);
Wwidth = sdpvar(N,1);
Wcenter = sdpvar(N,1);


x = sdpvar(5,1);

Q_mean_mat_var = sdpvar(4,N);
sb = sdpvar(N,1);
C = [0 0 0 0 1];

Y = [];
xk = x;

Q_mean_mat2 = Q_mean_mat_var;
set_back2 = sb;

F = [];

for k = 1:N

    sprintf('Building model progress %d %',round(k/N*100))
    W = Wwidth(k).*(lambda1(:,k)+lambda2(:,k));

    xk = Ad*xk - Bd(:,1)*V(k)+ Bd(:,2:end-1)*Q_mean_mat2(:,k) + Bd(:,end)*Wcenter(k) + W;

    Y = [Y;C*xk];
   
    F = [F, lambda1(:,k)-lambda2(:,k) == -(L(k,:)*kron(eye(1,N),Bd(:,1))')' + Bd(:,end), ...
        lambda1(:,k) >= 0, lambda2(:,k) >= 0];
   
    F = [F, V(k) +  (lambda3(k,:)+lambda4(k,:))*Wwidth <= 100, ...
        lambda3(k,:)-lambda4(k,:) == L(k,:), lambda3(k,:)>=0, lambda4(k,:) >= 0];
   
    F = [F, -V(k) +  (lambda5(k,:)+lambda6(k,:))*Wwidth <= 0, ...
        lambda5(k,:)-lambda6(k,:) == L(k,:), lambda5(k,:)>=0, lambda6(k,:) >= 0];
end


F = [F, 70 - slack - set_back2 <= Y <= 74 + slack + set_back2,  slack >= 0];%, 0 <= U <= 100];

objective = 10000000*sum(slack) + sum(V);

ops = sdpsettings('verbose',1, 'solver', '+cplex');
xk = [72.896;75.9814;75.8665;72.4428;72.5114];

set_back2 = set_back;
Q_mean_mat2 = Q_mean_mat;
Ureal = [];

controller = optimizer(F, objective,ops,{x(:,1), Q_mean_mat_var,sb,Wwidth,Wcenter},{V,L});

figure;
hold on;

for i = 1:48

    wn = rand(N,1);
    wc = (Q_lb(5,i:i+N-1)'+Q_ub(5,i:i+N-1)')/2;
    ww = (Q_ub(5,i:i+N-1)'-Q_lb(5,i:i+N-1)')/2;
    sol = controller{{xk(:,end), Q_mean_mat2(:,i:i+N-1),set_back2(i:i+N-1)',ww,wc}};
    L = sol{2};
    U = sol{1} + sol{2}*(Q_lb(5,i:i+N-1)' + wn.*(Q_ub(5,i:i+N-1)'-Q_lb(5,i:i+N-1)'));
%     uk = value(U(1));
    uk = U(1);
     xk = [xk Ad*xk(:,end) -  Bd(:,1)*uk + Bd(:,2:end-1)*Q_mean_mat(:,i) + Bd(:,end)*Q_real(5,i)];
    Ureal = [Ureal;uk];
   
    % Plot some random predictions
    %     if 0
    Y = [];
    clf;
   
    stairs(1:N,74 + set_back2(i:i+N-1)','LineWidth',4);hold on
    stairs(1:N,70 - set_back2(i:i+N-1)','LineWidth',4);
    for j = 1:100
        xp = xk(:,end-1);
        wn = rand(N,1);
        w = Q_lb(5,i:i+N-1)' + wn.*(Q_ub(5,i:i+N-1)'-Q_lb(5,i:i+N-1)');
        u = L*wn + V;

        for k = 1:N

            xp =  [xp Ad*xp(:,end) -  Bd(:,1)*U(k) + Bd(:,2:end-1)*Q_mean_mat(:,i+k-1) + Bd(:,end)*w(k)];
        end
        stairs(C*xp(:,2:end)); hold on
    end
    drawnow
    %     end
end

[Johan Löfberg]

I've tried to generalize the code slightly, so I think your case should work now. The robustification will run pretty slow as it is a quick hack to support this nonlinear paramterized case. Simply replace modules/robust/filter_normball

[Johan Löfberg]

Minor bug, use this

[LP]

Thank you for replying.


Yes it is working, but extremely (impossible for conducting simulation tests) slow.

What I would have tried now is to calculate the max part of min-max by hand (using strong duality), and call yalmip for only min part.

The approach is taken from this paper.

Rough description of the system is shown below.

For each constraint to be taken row wise (k), with bounds Omega, decision variable h and M, we have to perform:

max(||v|| <= Omega) F_k(A*x_0 + H*h + L*M*v + E*w) - f_k <= 0

max(||v|| <= Omega) S*M*v + S*h <= s

The only problem I am experiencing in this formulation (as you will probably see from the code too) is how to identify S and F matrix, since they are used in some constraints after taking the Kronecker product with an identity matrix(size of that identity is also not clear to me).

Due to the above mentioned problem, the M matrix shows NaN from second column onwards. And I am guessing because it never enters the optimization problem as a decision variable.

If you could kindly look into it, it would be great.

I will inform also, if I am able to solve it on my own.

Once again, thanks a lot for the help so far. Highly appreciated !

Regards,

[Johan Löfberg]

Applying full duality based robustification here is a waste and unnecessarily complicated (YALMIP would only do that if the uncertainty set was a general polytopic/conic set). For simple bounds, use that

f+g'*w>= 0 for all |w|<=r

is

f-r*norm(g,1)>= 0

(This is called the explicit filter in the YALMIP robustification module)

[Johan Löfberg]

Please supply the code which runs slow

[LP]

Thanks, attached "experiment_slow" is the main file which is run.

[Johan Löfberg]

This model appears to confuse YALMIP completely. As it is now, it is a model with high-level operators (norm) with  bilinearly parameterized arguments, and nonconvex use of the nonlinear operator (you have -norm), which effectively means YALMIP should introduce loads of binary variables. However, I only see 2 of them in the final model. 

Additionally, you have sums/differences of multiple nonlinear operators, all incluing uncertainty, which is really hard to robustify. If you would have had simple sums, YALMIP has measures to deal with this more or less conservatively (the auxreduce options), but in this nonconvex case, I have no idea what actually happens now.

It is not terribly slow though, takes 30 secs on my machine to derive the model. Considering the complexity of deriving (the probably wrong) model and the fact that it is a quick hack, I'm ok with that

[Johan Löfberg]

and now I see why YALMIP only introduces 2 binaries.

This

    objective = norm(1e6*slack(k),inf) + norm(U(k)*usep_in(k)*1800*0.00055,inf) - norm(R(k)*res_in(k)*1800*0.00055,inf);

    

is probably supposed to be

    objective = objective +norm(1e6*slack(k),inf) + norm(U(k)*usep_in(k)*1800*0.00055,inf) - norm(R(k)*res_in(k)*1800*0.00055,inf);

    

[LP]

Thanks a lot for the input. I will look into it for sure.

I think, that's why calculating max part by hand (explicitly) would be the best alternative, and call yalmip for minimizing

So I have changed the below equations

max(||v|| <= Omega) F_k(A*x_0 + H*h + L*M*v + E*w) - f_k <= 0

max(||v|| <= Omega) S*M*v + S*h <= s

to


F_k(A*x_0 + H*h)  - Omega*(norm(F_k*(L*Mk + Ek),1)) - f_k <= 0

S*h - Omega*(norm(S*M,1))  - s<= 0

But I get now the error saying

%%%%%%%%%%%%%%%ERROR%%%%%%%%%%%%%%%%%%%%%%
Error using sdpvar/model (line 63)
Failed when trying to create a model for the "norm" operator

Error in expandrecursive (line 19)
    [properties,F_graph,arguments,fcn] = model(variable,method,options,allExtStruct(ext_index),w);

Error in expandmodel>expand (line 376)
            [F_expand,failure,cause] =
            expandrecursive(recover(variables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'convex',allExtStruct,w);
           
Error in expandmodel (line 242)
                [F_expand,failure,cause] =
                expand(index_in_extended,variables,-sdpvar(F(constraint)),F_expand,extendedvariables,monomtable,variabletype,['constraint
                #' num2str(constraint)],0,options,
Error in compileinterfacedata (line 116)
    [F,failure,cause,operators] = expandmodel(F,h,options);

Error in export (line 135)
    [interfacedata,recoverdata,solver,diagnostic,F] = compileinterfacedata(F,[],logdetStruct,h,options,findallsolvers);

Error in optimizer (line 200)
    [aux1,aux2,aux3,model] = export((x == repmat(pi,nIn*mIn,1))+Constraints,Objective,options,[],[],0);

Files attached if you want to check it out what's going on real quick...

[Johan Löfberg]

Let me repeat my comment: Are you aware that the model you have is nonconvex (also if you remove all uncertainty and parameters)?

Why is it better to compute max manually (which is precisely what YALMIP does when it robustifies)?

I strongly doubt that your model is correct, as you don't appear to have any logic to deal with the complication that you have a sum of multiple norms involving uncertainty. It requires enumeration or projection, you cannot apply robustification on each term individually

[Johan Löfberg]

It doesn't even run here as the supplied prameters have the wrong dimension

>> {x1(:,1),x2(:,1), Q_mean_mat_var,sb,usep_in,res_in}

ans = 

  Columns 1 through 4

    [5x1 sdpvar]    [5x1 sdpvar]    [4x48 sdpvar]    [48x1 sdpvar]

  Columns 5 through 6

    [48x1 sdpvar]    [48x1 sdpvar]

>> {xk1(:,end), xk2(:,end), Q_mean_mat2(:,i:i+N-1),set_back2(i:i+N-1)',usep(i:i+N-1)',res(i:i+N-1)'}

ans = 

  Columns 1 through 4

    [5x1 double]    [5x1 double]    [4x48 double]    [48x1 double]

  Columns 5 through 6

    [1x48 double]    [1x48 double]

[LP]

Ok, this I don't understand.

Let's say for the deterministic case, removing all the uncertainties and approximation. The most basic system is represented as:

min  (U) + min(slack)

s.t.
 
x = Ax + Bu + Ew
Cx - AA - slack <= 0
Cx +AA + slack <= 0
slack >= 0
0 <= u <= 100

The terms included afterwards e.g. (price, revenues), though time-varying but are assumed to be known.

[Johan Löfberg]

..and you appear to be working with two completely different models when you send codes. Very confusing, so my comments only apply to the YALMIP model with nonconve nobjective. The model where you try to do things manually, I am not looking at at all, as there is no chance that I would be able to see what you are doing there, 

[LP]


Sorry for the mistake, the transpose of these two terms should be removed. --> usep(i:i+N-1)',res(i:i+N-1)'

[Johan Löfberg]

Your model has ...+norm(U,inf)-norm(R,inf), hence nonconvex and requires MILP modelling

[Johan Löfberg]

This can be fixed by adding

X = reshape(X,[],1);

Z = reshape(Z,[],1);

on line 107 in sdpvar/norm

Thank you for alerting me

[LP]

Thank you.

Have you checked solving also, after the model has been created ?

for me it is taking a lot of time with cplex ... 

any idea of different solver who might be a bit faster ?

[Johan Löfberg]

Which model are you talking about. As I said, you have two conceptually very different models, one convex and one nonconvex, and both of them had errors making them impossible to run, and one of them having terms missing in the objective..

[LP]

Ok about the two models. They are getting nonconvex and convex due to the operators being used by me naively. And possible because difference lies in the objective function.
One only minimizes U and slack. The other minimizing the cost and maximizes the reserves so it has U.*usep and R.*res

the system matrices (A,B,C, E) for both are the same.

In the one minimizing U and slack, I try to augment the matrices and do the following

max(||v|| <= Omega) F_k(A*x_0 + H*h + L*M*v + E*w) - f_k <= 0

max(||v|| <= Omega) S*M*v + S*h <= s

The other minimizing the cost and maximizes the reserves so it has U.*usep and R.*res I try to use "robustify"

[Johan Löfberg]

so which model are you asking about

[LP]

the one with the robustify. I have attached it. please let me know if it is atleast compiling.

Sorry for the confusion, and thanks a million for being patient.

Regards,

[Johan Löfberg]

Sure, it creates the controller in 30s or so on my machine. You get a horrible model though which takes ages to solve. Not much you can do there, you've formulated a bad model when you started adding nonconvex terms to the objective