infeasible?

[Ali Esmaeilpour]

Hello professor
again I'm facing not feasible
here is code:

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;

Ek1 = zeros(N+1,1);
Ek2 = ones(N+1,1);
Ek3 = ones(N,1);
Ek4 = zeros(N,1);

constraint=[];
 
objective=0;

for k = 1:N
   
    if k==1
        h(:,:) = [kron(Ek2,hx);kron(Ek3,hu)];
    else
        h(:,:) = [kron(Ek1,hx);kron(Ek4,hu)];
    end
   
   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*(sqrt(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h))<=g, (h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h)>=0, (h'*(Abar+Bhat*K{k}*Ek))>=0,((Abar+Bhat*K{k}*Ek)'*h)>=0, (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')>=0;
end

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

K = value(K);

[Johan Löfberg]

First, your code does not run as it is missing a closing bracket on the constraint

Second, it is trivially infeasible already when you construct it, as (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw') evaluates to a fixed matrix, and then you cannot write ...>=0 as it simply will be a large matrix of true/false

[Ali Esmaeilpour]

ok tnx but i wanna ask if those constraints are in a matrix form should i write it exactly as it is or i can put them one by one i mean get them out of matrix and write as >=0 for each

constraints are like the picture i attached

[Ali Esmaeilpour]

doing those things i said i still got error:
"One of the constraints evaluates to a FALSE LOGICAL variable"

is it because of this??? (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')

code is as follows:

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

[Johan Löfberg]

When N is 5, 

[(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h) (h'*(Abar+Bhat*K{k}*Ek));((Abar+Bhat*K{k}*Ek)'*h) (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')]

is a fixed matrix (my previous statement was due to the strange indentation / missing ] in your code). Hence, it makes no sense to try to add it to the model

>> [(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h) (h'*(Abar+Bhat*K{k}*Ek));((Abar+Bhat*K{k}*Ek)'*h) (Gxx*Sigma0*Gxx' + Gxw*Sigmaw*Gxw')]

ans =

         0         0         0         0         0         0         0         0         0         0         0

         0    2.0206    0.0000    2.0413    0.0020    2.0622    0.0041    2.0834    0.0062    2.1047    0.0084

         0    0.0000    2.0198   -0.0020    2.0397   -0.0041    2.0598   -0.0062    2.0800   -0.0084    2.1005

         0    2.0413   -0.0020    3.0623    0.0000    3.0936    0.0031    3.1254    0.0063    3.1574    0.0095

         0    0.0020    2.0397    0.0000    3.0598   -0.0031    3.0899   -0.0062    3.1203   -0.0095    3.1511

         0    2.0622   -0.0041    3.0936   -0.0031    4.1254    0.0000    4.1677    0.0042    4.2104    0.0084

         0    0.0041    2.0598    0.0031    3.0899    0.0000    4.1203   -0.0042    4.1609   -0.0084    4.2019

         0    2.0834   -0.0062    3.1254   -0.0062    4.1677   -0.0042    5.2104    0.0000    5.2638    0.0053

         0    0.0062    2.0800    0.0063    3.1203    0.0042    4.1609    0.0000    5.2019   -0.0052    5.2531

         0    2.1047   -0.0084    3.1574   -0.0095    4.2104   -0.0084    5.2638   -0.0052    6.3178    0.0000

         0    0.0084    2.1005    0.0095    3.1511    0.0084    4.2019    0.0053    5.2531    0.0000    6.3049

same thing with the other expression, all terms are fixed, so you simply evaluate a comparison

h'*(Cs+Bhat*K{k}*Euhat)+r*(sqrt(h'*(Abar+Bhat*K{k}*Ek)*M*(Abar+Bhat*K{k}*Ek)'*h))<=g

ans =

  1×20 logical array

   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1

[Ali Esmaeilpour]

oh there is K{k} as a decision variable and btw it's there still we got logical values?

[Johan Löfberg]

so what? your data causes your prediction to be constant, such as

>> ((Abar+Bhat*K{k}*Ek)'*h)

[Ali Esmaeilpour]

I think I found the cause of getting logical variables i've put an if else there which depends on k=0,..,N and i=0,..,m and its form is in the picture. am I right? and if that is the problem how can i put it in my loop to have decision variables involved instead of logical ones???

[Ali Esmaeilpour]

that h has to depend on k but it is zero and cause this to be zero. should i define it as a sdpvar?

[Johan Löfberg]

No. You simply have to understand why it happens to be that h is such that ((Abar+Bhat*K{k}*Ek)'*h) is 0. If it is intentional, then it makes no sense to add any constraint for k=N, as everything is fixed anyway. If it is not intentional, well then you have to sort out the bug you have in defining the data

[Ali Esmaeilpour]

ok i've added this:

  if k==1
        h(:,:) = [kron(Ek2,hx);kron(Ek3,hu)];
    else
        h(:,:) = [kron(Ek1,hx);kron(Ek4,hu)];
    end

and i mean for k=1 h is a different matrix and it is not zero but after running my code h is zero and i'm not sure if it's a different matrix for k=1
how can i figure it out?

[Johan Löfberg]

That's not yalmip related but very specific to your problem, and something you will have to figure out.

[Ali Esmaeilpour]

ok, tnx professor for your answers