Model creation failed

[Ali Esmaeilpour]

Hello professor Lofberg
I've problem wit the following code. It says "Model creation failed" and see https://yalmip.github.io/squareroots/ but after reading that article I couldn't figure out how to debug my code.
my code is:

clc;
clear;
close all;
%% Time structure
dt=0.01;          %sampling time

%% System structure
A = [1.02  -0.1
     0.1   0.98];
 
B = [0.5 0
    0.05 0.5];

G = [0.3 0
     0 0.3];
 
Qf = [50 0
      0 50];

[A, B] = c2d(A,B,dt);

[nx,nu] = size(B);

%% MPC parameters
Q  = eye(2);

R  = eye(2)*50;

N = 5;

%% Building block matrices
I = eye(2,2);

q = 2;

m = 2;

Qbar = blkdiag(Q,Q,Q,Q,Qf,R,R,R,R);

Fq = blkdiag(sqrt(Q),sqrt(Q),sqrt(Q),sqrt(Q),sqrt(Q),sqrt(Qf),sqrt(R),sqrt(R),sqrt(R),sqrt(R),sqrt(R));

Gxx = zeros(2*N , size(A,1));
for i = 1:N
    Gxx(q*i-q+1:q*i , :) = A^i;
end

Gxu = zeros(2*N , 2*N);
for i = 1:N
   for j = 1:i
      Gxu(q*i-q+1:q*i , m*j-m+1:m*j) = A^(i-j)*(B);
   end   
end

Gxw = zeros(2*N , 2*N);
for i = 1:N
    for j = 1:i
    Gxw(q*i-q+1:q*i , m*j-m+1:m*j) = A^(i-j)*(I);   
    end
end

%% Solve using YALMIP

K = sdpvar(repmat(10,1,N),repmat(22,1,N));

alpha = 0.975;

r = 1.96;

Sigmaw = eye(10,10);

Sigma0 = eye(2,2);

xbar0 = ones(2,20);

Cs = [Gxx*xbar0;zeros(12,20)];

hx = [-1/sqrt(5);-2/sqrt(5)];

hu = [0;0];

Ahat = [Gxx*xbar0 eye(10,2);zeros(12,20) zeros(12,2)];

Bhat = [Gxu;eye(12,10)];

Abar = [eye(11,10);zeros(11,10)];

M = Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw';

Mhat = [M zeros(10,12);zeros(12,10) ones(12,12)];

Euhat = ones(22,20);

Ek = [zeros(12,10);eye(10,10)];

g = 3;

constraint=[];
 
objective=0;

h = [-0.447213595499958;-0.894427190999916;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];

for k = 1:N  

b = sqrt(h'*((Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)')*h);

objective = objective + 0.5*trace(Fq'*(Ahat+Bhat*K{k})*Mhat*(Ahat+Bhat*K{k})'*Fq);

constraint = [constraint, ((h'*(Cs+Bhat*K{k}*Euhat))+(r*b))<=g, [b^2 (h'*(Abar+Bhat*K{k}*Ek));((Abar+Bhat*K{k}*Ek)'*h) (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')]>=0];

end

options = sdpsettings('solver','sedumi');
optimize (constraint,objective,options);

K = value(K);

[Johan Löfberg]

sqrt is a convexity aware operator, and will fail as sqrt(quadratic) does not satisfy rules for convexity propagation

You should simply write it as norm((chol(M)*(Abar+Bhat*K{k}*Ek)')*h)

However, it will still fail as your semidefinite constraint involving b^2 won't be linear or convex, and you thus have no solver to solve it. 

[Ali Esmaeilpour]

ok, but "b" showed up as is in the picture I attached. b=sqrt(h'*sigmaz*h)
I dont know that I'm wrong about recognizing "b" or not?

[Johan Löfberg]

I have no idea what you are saying. Just because some paper write an expression in some form, does not mean it is the form that can be used in modelling language

the derivation is wrong, as there should be an inverse on M in the final expression

in the derivation, b is simply a new auxiliary variable, intended to upper bound the norm. nevertheless, the final expression obtained is not an LMI, as it involves b^2

[Ali Esmaeilpour]

eh... excuse me there should be inv(M)?

another subject is that you mean I can't say b=sqrt(h'*sigmaz*h)?

[Johan Löfberg]

b is a new variable, intended to act as an upper bound on sqrt(h'*sigmaz*h). EQuality will only hold at optimality, if required

[Ali Esmaeilpour]

so there is no way to figure out how to use b in modeling language if we assume optimality?

[Johan Löfberg]

??

I already showed you how to model the original stuff. Skip the whole failed Schur complement stuff and reformulations and stuff. You simply have h_k^T*Z + r*norm(chol(M)*(A+B*K)*h_k) <= g