Re: there is a problem in example of LPV-MPC in yalmip

[Johan Löfberg]

obj is defined during the "First step" part of the code. You have to copy-paste various blocks to setup the complete logic

[asemeh k]

thanks Prof for your answer

sorry for bothering you,but I have another question about paper "Explicit Model Predictive Control for Linear Parameter-Varying Systems".

in numerical example how can obtain A1 ,A2,B1, B2  and scheduling parameter from this system.

best

Asi

[Johan Löfberg]

I guess some outer bounding of nonlinearity

[asemeh k]

Dear Prof Johan 

Thank you again for your answering .

when run the code (all steps in on code ) I find  problem again about intermediate step and final step like as 

Index exceeds matrix dimensions.

Error in yalmip1 (line 104)

 S = unique(([reshape([sol{k+1}{1}.Bi{:}]' ,nx,[])'  reshape([sol{k+1}{1}.Ci{:}]'

 ,[],nu)]),'rows');

and if I run code include just intermediate step I find this error  

Undefined variable "sol" or class "sol".

Error in yalmip_intermediat (line 61)

 S = unique(([reshape([sol{k+1}{1}.Bi{:}]' ,nx,[])'  reshape([sol{k+1}{1}.Ci{:}]'

 ,[],nu)]),'rows');

 

I have another questions, 

How can I obtain figure of result from  paper ? I mean which code I should use for plot this figure.

Best

Asi

[Kolapo Alli]

Is there a way of representing this convolution expression below in matrix notation?

int(z(y)*(int(f_w(y)*g_w(y)'*(z(y)-w(y)),0,wr)),0,wr)

 where f_w=sdpvar(n,1); g_w=sdp(n,1). The above solution expression actually ran well but I have a concern on how to handle the integral parts which I think should be handle in matrix format, since problems in sdp are written in matrix?

[Johan Löfberg]

I have no idea what you are trying to ask, or do and the code doesn't make sense. If f_w is a variable, what do you mean with f_w(y). Is y a number between 1 and n?

[Kolapo Alli]

Ofcourse y=1:n

int(z(y)*(int(f_w(y)*g_w(y)'*(z(y)-w(y)),0,wr)),0,wr)

This expression was actually a part of the problems that I am solving in my research work using sdp technique. I want to know how the integral parts can be handle in matrix notation.

[Johan Löfberg]

Still code doesn't make sense. The second argument to int should be a variable indicating the variable which is integrated, in your code you have a 0

[Kolapo Alli]

I actually got good result from the code but my supervisor insisted that there should be a way of reformulating the integrals in sdp formation so as to show to every readers that the problem was actually be reformulated in matrix notation.

[Johan Löfberg]

You would have to supply code which can be run to actually show what you are doing

[Kolapo Alli]

Sooner, I will get back to you on the complete code that I am working on.  

[Ros Cibely]

Hello, I'm traying reproduce the example of YALMIP site: https://yalmip.github.io/example/explicitlpvampc/, and I would like to plot the response of the system.
The variable x which has the model states is NaN. Why?

This is the code:

close all, clear all, clc

yalmip('clear')

% Model data
A1 = [ 2.1 1;0.3 0];
A2 = [-2.1 1;0.3 0];
B = [1;0.1];
f = [1;0];

% System sizes
nx = 2; % Number of states
nu = 1; % Number of inputs
ndyn = 2; % Number of vertex systems

% State and input constraints
xmin = [-10;-10];
xmax = [ 10; 10];
umin = -5;
umax = 5;

% Reference point
xref = [0.63;0.19];

% MPC data
Q = eye(nx);
R = 0.1;
N = 4;
nrm = 1;

% States x(k), ..., x(k+N)
x = sdpvar(repmat(nx,1,N),repmat(1,1,N));
% Inputs u(k), ..., u(k+N) (last one not used)
u = sdpvar(repmat(nu*ndyn,1,N),repmat(1,1,N));
% Scheduling parameter
th = binvar(repmat(ndyn,1,N),repmat(1,1,N));
% Epigraph variable
sdpvar w

for k = N:-1:1 % shifted: N-1:-1:0

% Parameter simplex
F = [uncertain(th{k}), sum(th{k}) == 1, 0 <= th{k} <= 1];
F = [F, 0 <= w <= 10000];

% Uncertain predictions and control
uth = kron(th{k},eye(nu))'*u{k}; % u(th)
Ath = [A1 A2]*kron(th{k},eye(nx)); % A(th) = A1*th1 + A2*th2 + ...
xp = Ath*x{k} + f + B*uth;

% Input constraints
F = [F, repmat(umin,ndyn,1) <= u{k} <= repmat(umax,ndyn,1)];

% State constraints
F = [F, xmin <= x{k} <= xmax];

% Insert step-specific code here
% First-step
% |x| written as max(a[x;u]+b)
a = [-1+2*dec2decbin(0:2^nx-1,nx) zeros(2^nx,nu)];
b = zeros(2*nx,1);
% |u| written as max(c[x;u]+d)
c = [zeros(2^nu,nx) -1+2*dec2decbin(0:2^nu-1,nu)];
d = zeros(2*nu,1);

% |x|+|u| written as max(aa[x;u]+bb)
aa = repmat(a,2*nu,1) + kron(c,ones(size(a,1),1));
bb = repmat(b,2*nu,1) + kron(d,ones(size(a,1),1));

% Final state constraints
F = [F, xmin <= xp <= xmax];

% Cost function
F = [F, aa*[Q*(xp-xref);R*uth]+bb + norm(Q*(x{k}-xref),nrm) <= w];
obj = w;

% Determine robust counterpart
[F,obj] = robustify(F,obj,[],th{k});

% % Intermediate Step
% % Get the hyperplanes for cost-to-go
% S = unique(([reshape([sol{k+1}{1}.Bi{:}]' ,nx,[])' reshape([sol{k+1}{1}.Ci{:}]' ,[],nu)]),'rows');
% a = [S(:,1:nx) zeros(size(S,1),nu)];
% b = S(:,nx+1:end);
%
% % Jp+|u| = max(a[x;u]+b)+|u|=max(a[x;u]+b)+max(c[x;u]+d)=max(aa[x;u]+bb)
% aa = repmat(a,2*nu,1) + kron(c,ones(size(a,1),1));
% bb = repmat(b,2*nu,1) + kron(d,ones(size(a,1),1));
%
% % Constrain predicted state
% [H,K] = value(sol{k+1}{1}.Pfinal);
% F = [F, H*xp <= K];
%
% % Cost function
% F = [F, aa*[xp;R*uth]+bb + norm(Q*(x{k}-xref),nrm)< w];
% obj = w;
%
% % Determine robust counterpart
% [F,obj] = robustify(F,obj,[],th{k});

% % Final Step
%
% % Get the hyperplanes for cost-to-go
% S = unique(([reshape([sol{k+1}{1}.Bi{:}]',nx,[])' reshape([sol{k+1}{1}.Ci{:}]',[],nu)]),'rows');
% a = [S(:,1:nx) zeros(size(S,1),nu)];
% b = S(:,nx+1:end);
%
% % Jp+|u| = max(a[x;u]+b)+|u|=max(a[x;u]+b)+max(c[x;u]+d)=max(aa[x;u]+bb)
% aa = repmat(a,2*nu,1) + kron(c,ones(size(a,1),1));
% bb = repmat(b,2*nu,1) + kron(d,ones(size(a,1),1));
%
% % Redefine final step input
% u{1} = sdpvar(nu,1); % final step input
%
% % Input constraints
% F = [umin <= u{1} <= umax];
%
% % State update equation
% xp = x{k} + f + B*u{1}; % x{1} = z
%
% % Constrain predicted state
% [H,K] = value(sol{k+1}{1}.Pfinal);
% F = [F, H*xp <= K];
%
% % Cost function
% F = [F, aa*[xp;R*u{1}]+bb <= w];
% obj = w;
%
% Solve multi-parametric problem
 [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,[],x{k},u{k});
end

The variável sol has the solution, hasn't?

Could someone help me with this doubts?

[Johan Löfberg]

so how have you used the computed solution to perform a simulation to create the states? The code you supply only shows that you generate the solution, and that code appears to work as a solution map was gerenated

>> assign(x{1},[0.1;0.2]);

>> value(Optimizer{1})

ans =

   -0.7800

   -0.4300

Yes, sol contains the PWA solution in low-level MPT format

[Ros Cibely]

Thanks! So, I have to use the cell sol to perform a simulation to create the response states?

[Johan Löfberg]

or use last output Optimizer as I did

[Ros Cibely]

Where did you find the vector value [0.1;0.2]? You use it in assign(x{1},[0.1;0.2])?


[Johan Löfberg]

just an example

[Ros Cibely]

How I can use the cell sol to perform a simulation to create the states?

[Ros Cibely]

I know that the region of Pn containing the given state x(k), but sol{1}{1}.Pn is a polyhedral

[Johan Löfberg]

huh?

[Ros Cibely]

Well, I would like to know how to plot the output of the system model given..

[Johan Löfberg]

>> [~,where]=isinside(sol{1}{1}.Pn,x0);

>> sol{1}{1}.Fi{where}*x0+sol{1}{1}.Gi{where}

ans =

   -0.7800

   -0.4300

[Ros Cibely]

Thank you very much Johan Löfberg.

[Ros Cibely]

Hello Johan,

I am trying to compile all the code with the frist step and the intermediete step, but I think there is a error in the code example

Model predictive control - LPV models (https://yalmip.github.io/example/explicitlpvmpc/)

Error using value
Too many output arguments.

Error in MPC_LPV (line 112)

    [H,K] = value(sol{k+1}{1}.Pfinal);

value(sol{k+1}{1}.Pfinal) I believe this command should return numeric values, but it doesn't.


%%
%https://yalmip.github.io/example/explicitlpvmpc/
%%
clc, clear
% YALMIP options
yalmip('clear')
yopts = sdpsettings('robust.polya',1);

% Model data
A1 = [0.85 0;0.25 0.65];
A2 = [0.85 0;-0.3 0.65];
B1 = [1;-1];
B2 = [1;1];
B{1} = B1;
B{2} = B2;

% System sizes
nx   = 2; % Number of states
nu   = 1; % Number of inputs
ndyn = 2; % Number of vertex systems

% State and input constraints
xmin = [-10;-10];

xmax = [  8;  8];
umin = -0.5;
umax =   1 ;

% MPC data
Q   = eye(nx);

R   = 0.01;
N   = 3;
nrm = inf;
xref = [0;0];

% States x(k), ..., x(k+N)
x = sdpvar(repmat(nx,1,N),repmat(1,1,N));
% Inputs u(k), ..., u(k+N) (last one not used)
u = sdpvar(repmat(nu*ndyn,1,N),repmat(1,1,N));
% Scheduling parameter
th = binvar(repmat(ndyn,1,N),repmat(1,1,N));
% Epigraph variable
sdpvar w;

for k = N:-1:1   % shifted: N-1:-1:0
    % Parameter simplex
    F = [uncertain(th{k}), sum(th{k}) == 1, 0 <= th{k} <= 1];
    F = [F, 0 <= w <= 10000];
    % Uncertain predictions and control
    uth = kron(th{k},eye(nu))'*u{k};    % u(th)
    Ath = [A1 A2]*kron(th{k},eye(nx));  % A(th) = A1*th1 + A2*th2 + ...

    Bth = [B1 B2]*kron(th{k},eye(nu));  % B(th)
    xp  = Ath*x{k} + Bth*uth;

    % Input constraints
    F = [F, repmat(umin,ndyn,1) <= u{k} <=  repmat(umax,ndyn,1)];
    % State constraints
    F = [F, xmin  <= x{k} <= xmax];
   
    % Insert step-specific code here

    %%
    %First step

    % |x| written as max(a[x;u]+b)

    a = [kron(eye(nx),[1 -1]') zeros(2*nx,nu)];

    b = zeros(2*nx,1);
    % |u| written as max(c[x;u]+d)

    c = [zeros(2*nu,nx) kron(eye(nu),[1 -1]')];

    d = zeros(2*nu,1);
    % |x|+|u| written as max(aa[x;u]+bb)
    aa = repmat(a,2*nu,1) + kron(c,ones(size(a,1),1));
    bb = repmat(b,2*nu,1) + kron(d,ones(size(a,1),1));
   
    % Final state constraints
    F = [F, xmin  <= xp   <= xmax];
    % Cost function
    F   = [F, aa*[Q*(xp-xref);R*uth]+bb + norm(Q*(x{k}-xref),nrm) <= w];
    obj = w;
    % Determine robust counterpart

    [F,obj] = robustify(F,obj,yopts,th{k});
    %%
   
   
    % Solve multi-parametric problem
    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,yopts,x{k},u{k});
end

for k = N-1:-1:1   % shifted: N-1:-1:0

    % Parameter simplex
    F = [uncertain(th{k}), sum(th{k}) == 1, 0 <= th{k} <= 1];
    F = [F, 0 <= w <= 10000];
    % Uncertain predictions and control
    uth = kron(th{k},eye(nu))'*u{k};    % u(th)
    Ath = [A1 A2]*kron(th{k},eye(nx));  % A(th) = A1*th1 + A2*th2 + ...

    Bth = [B1 B2]*kron(th{k},eye(nu));  % B(th)
    xp  = Ath*x{k} + Bth*uth;

    % Input constraints
    F = [F, repmat(umin,ndyn,1) <= u{k} <=  repmat(umax,ndyn,1)];
    % State constraints
    F = [F, xmin  <= x{k} <= xmax];
   
    % Insert step-specific code here

    %%
    %Intermediate steps

    % Get the hyperplanes for cost-to-go

    S = unique(([reshape([sol{k+1}{1}.Bi{:}]' ,nx,[])'  reshape([sol{k+1}{1}.Ci{:}]' ,[],nu)]),'rows');

    a = [S(:,1:nx) zeros(size(S,1),nu)];

    b = S(:,nx+1:end);
   

    % Jp+|u| = max(a[x;u]+b)+|u|=max(a[x;u]+b)+max(c[x;u]+d)=max(aa[x;u]+bb)

    aa = repmat(a,2*nu,1) + kron(c,ones(size(a,1),1));
    bb = repmat(b,2*nu,1) + kron(d,ones(size(a,1),1));
   

    % Constrain predicted state

    [H,K] = value(sol{k+1}{1}.Pfinal);

    F     = [F, H*xp <= K];
   

    % Cost function
    F   = [F, aa*[xp;R*uth]+bb <= w];
    obj = norm(Q*(x{k}-xref),nrm) + w;
   
    %%
   
   
    % Solve multi-parametric problem
    [sol{k},diagnost{k},Uz{k},J{k},Optimizer{k}] = solvemp(F,obj,yopts,x{k},u{k});
end

[Ros Cibely]

Can I use double in place value?
I mean double(sol{k+1}{1}.Pfinal)

[Johan Löfberg]

I think that is the command, yes (Pn is not a YALMIP object, it is a polytopic object from MPT)

[Johan Löfberg]

your code is weird incorrect though, as you perform the initial first step N times, the intermediate step N times, and then never perform the final step.

[misagh ahmad]

Hello
excuse me,
I run following code for optimize a H controller, but the result is "NAN", while the LMI structure in my problem have no problem, I try two solve with both rolmip & yalmip but no answer except "NAN"
can you help me please
clear all
clc
 m=ureal('m',1.5,'Range',[1 2]);
l=ureal('l',1.2,'Range',[.8 1.4]);
R=.25 ,M=55+m , D=.15;
g=9.8;
Jsy=1.5,Jp=m*l^2  ,Jfi=2.5+m*l^2 ,omga=Jfi*M-.9*m^2*l^2 ;

taw=.5;
alfa=.25;
%syms x1 x2 x3 x4 x5 x6
%x=[x1 x2 x3 x4 x5 x6]'
 
A=[0 1 0 0 0 0;0 0 -.9*((m*l)^2*g)/omga 0 0 0;0 0 0 1 0 0;0 0 -M*m*g*l/omga 0 0 0;0 0 0 0 0 1;0 0 0 0 0 0];
A=A+eye(6)

B2=Jfi/R/omga*[0 0;1 1;0 0;omga omga;0 0;D*omga/2/Jsy/Jfi -D*omga/2/Jsy/Jfi];
B1=1/omga*[0 0 0 0;Jfi Jfi 0 0;0 0 0 0;-.95*m*l -.95*m*l 0 0;0 0 0 0;D*omga/2/Jsy -D*omga/2/Jsy 0 0];

C2=[0 1 0 0 0 0
    0 0 1 0 0 0
    0 0 0 1 0 0
    0 0 0 0 1 0];
C1=[0 0 1 0 0 0
    0 1 0 0 0 0
    0 0 0 1 0 0
    0 0 0 0 1 0];
D=[0 0 0 0
    0 0 -1 0
    0 0 0 0
    0 0 0 -1];
H=[0 0 0 0 0 0
    0 0 0 0 0 0
    0 0 1 1 0 0
    0 0 1 1 0 0
    0 0 0 0 0 0
    0 0 0 0 0 0];
taw=1;
alfa=1;
X=sdpvar(6,6);
Y=sdpvar(2,6);
gma=sdpvar(1,1);
for(i=1)
A=usample(A,1);
B1=usample(B1,1);
B2=usample(B2,1);
pscy=[X*A'+Y'*B2'+A*X+B2*Y B1 eye(6,6) X*C1' X*H';
B1' -gma^2.*eye(4,4) zeros(4,6) zeros(4,4) zeros(4,6);
eye(6,6) zeros(6,4) -taw.*eye(6,6) zeros(6,4) zeros(6,6);
C1*X zeros(4,4) zeros(4,6) -eye(4,4) zeros(4,6);
H*X zeros(6,4) zeros(6,6) zeros(6,4) (-1/taw/alfa^2).*eye(6,6)];
end
constraint = [X>=0, pscy<=0];
    objective =gma;
 
    sol = optimize(constraint,objective)
X=value(X)
Y=value(Y)
K=pinv(C)*pinv(Y);

[Johan Löfberg]

you cannot use ureal with YALMIP. That is something from the robust control toolbox or something like that

[misagh ahmad]

OK,that would be solved
 but how can I calculate the following problem,
X=sdpvar(6,6);
Y=sdpvar(2,6);

C= defined matrix with proper dimension

 where in main problem Y=KCX and my goal is to calculate variable of K=YX^-1C^-1, while there is no strict way to invert X

thanks

​

[Johan Löfberg]

If you have to work optimize over an inverse, and cannot come up with a linerizing convexifying variable change or congruence etc, you have to accept that it is a very (very) hard problem. Common trick is to replace X^-1 with a new variable Q, and add the nonconvex nonlinear equality X*Q == eye(n).  If you are doing standard state feedback, the are classical variable changes, congruences and Schur complement tricks

[舒彤]

Dear Prof Johan and fati-a:

 there is a problem in example of LPV-MPC in yalmip

***** Starting YALMIP robustification module. *********************

 

- Detected 2 uncertain variables

 - Detected 1 independent group(s) of uncertain variables

 - Enumerated 2 vertices

***** Derivation of robust counterpart done ***********************

Undefined variable "sol" or class "sol".

Error in lpvmpc (line 79)

S = unique(([reshape([sol{k+1}{1}.Bi{:}]' ,nx,[])'  reshape([sol{k+1}{1}.Ci{:}]' ,[],nu)]),'rows'); 

I have same questions, 

How can I solve this problem?

Best

shutong

[Johan Löfberg]

you would have to supply the full code you are running

[Johan Löfberg]

looking t the error message, it simply looks like you haven't constructed a correct version of the complete code, as it says that there is no variable sol available, and sol is generated when you compute the parametric solution so it cannot be missing unless you've missed copy-pasting some code