Constraints problem

[optimizer]

Hi Johan, 

I am using the function you can find below.

When i run the optimizer without ConstrS using the debug function i see this result about xU and xL constraints and others:

>> [model.lb model.x0 model.ub]

ans =

            0       1.0451         1.08

            0      0.38147         1.08

            0     0.053452         1.08

            0      0.15835         1.08

            0      0.80087         1.08

            0        75153    1.449e+05

            0   1.0215e+05    1.449e+05

            0        68801    1.449e+05

            0         4964    1.449e+05

            0        72362    1.449e+05

      -0.5455         -0.2      -0.1818

      -0.0364        -0.03      -0.0182

       -0.019        -0.01      -0.0019

      -0.9524         -0.6       -0.381

         -Inf     -0.11579          Inf

         -Inf   -0.0044277          Inf

         -Inf       -41465          Inf

         -Inf      -167.12          Inf

         -Inf   -0.0056313          Inf

         -Inf      -0.1006          Inf

         -Inf      -713.07          Inf

         -Inf      -3746.3          Inf

  -1.8484e+05            0        93979

         -Inf            0          Inf

       -16356            0        23361

         -Inf            0          Inf

      -5815.2            0       7145.6

      -1758.2            0       959.01

         -Inf            0          Inf

       -28340            0        32980

         -Inf            0          Inf

  -2.1169e+05            0   1.3436e+05

And everything is fine to me.

While when i run the optimizer with ConstrS as another constraint i got:

>> [model.lb model.x0 model.ub]

ans =

            0       1.0451      0.34207

        5e-16      0.38147        5e-16

        5e-16     0.053452        5e-16

        5e-16      0.15835        5e-16

        5e-16      0.80087        5e-16

            0        75153        34117

        5e-16   1.0215e+05        5e-16

        5e-16        68801        5e-16

        5e-16         4964        5e-16

        5e-16        72362        5e-16

      -0.5455         -0.2      -0.1818

      -0.0364        -0.03      -0.0182

       -0.019        -0.01      -0.0019

      -0.9524         -0.6       -0.381

         -Inf     -0.11579          Inf

         -Inf   -0.0044277          Inf

         -Inf       -41465          Inf

         -Inf      -167.12          Inf

         -Inf   -0.0056313          Inf

         -Inf      -0.1006          Inf

         -Inf      -713.07          Inf

         -Inf      -3746.3          Inf

  -1.8484e+05            0        93979

         -Inf            0          Inf

       -16356            0        23361

         -Inf            0          Inf

      -5815.2            0       7145.6

      -1758.2            0       959.01

         -Inf            0          Inf

       -28340            0        32980

         -Inf            0          Inf

  -2.1169e+05            0   1.3436e+05

The problem is that i don't want to modify the first 10 constraints while, as you can see xU and xL constraints are changed by introducing the matricial constraint ConstrS. Why? How do you suggest to handle it?

This is my function:

function [DA,DC] = my_function(A,C,DA,DC,U_s,o_s,x,xL,xU,o,y,yL,yU,U,xs)

n=size(A,1);

n_ep=length(o); 

p=cell(n_ep,1); 

for i=1:n_ep    

    p{i}=(o(i).^2*(A+DA)-(C+DC))*U(:,i); 

end

pvec=vertcat(p{:});    

B=getbase(pvec);        

v=getvariables(pvec);   

Q=eye(size(B,2));       

idx_x=depends(x).';     

idx_y=depends(y).';     

for ii=1:length(v)

    [~,~,i_x]=intersect(depends(recover(v(ii))),idx_x);

    if ~isempty(i_x)

        Q(ii+1,ii+1)=1/xs(i_x);       

    end

end

B=B*Q; 

pvec = B*[1;recover(v)];

v=getvariables(pvec);

v1=setdiff(getvariables(pvec),depends(pvec)); 

t=sdpvar(length(v1),1);                      

v2=getvariables(t);                     

numt=length(v2);                        

w=sdpvar(length(v),1);                 

LIA=ismember(v,v1);                     

v(LIA)=v2;                            

rvec = B*[1;recover(v)];                 

s = sdpvar(n*n_ep,1);    

xL=xL.*xs; 

xU=xU.*xs; 

ConstrT=[];         

ConstrTRelax=[];    

for i=1:numt        

    [~,~,i_x]=intersect(depends(recover(v1(i))),idx_x);

    [~,~,i_y]=intersect(depends(recover(v1(i))),idx_y);

    t_i=v2(i);

    ConstrT=ConstrT+[t(i)==x(i_x)*y(i_y)];

    ConstrTRelax=ConstrTRelax+[recover(t_i) >= xL(i_x)*y(i_y) + yL(i_y)*x(i_x) - xL(i_x)*yL(i_y)]+...

                              [recover(t_i) >= xU(i_x)*y(i_y) + yU(i_y)*x(i_x) - xU(i_x)*yU(i_y)]+...

                              [recover(t_i) <= xU(i_x)*y(i_y) + yL(i_y)*x(i_x) - xU(i_x)*yL(i_y)]+...

                              [recover(t_i) <= xL(i_x)*y(i_y) + yU(i_y)*x(i_x) - xL(i_x)*yU(i_y)];

end

ConstrX=[xL<=x<=xU]:'X Feasibility';

ConstrY=[yL<=y<=yU]:'Y Feasibility';

ConstrBounds=[ConstrX,ConstrY];   

ConstrS=[];

ConstrS=ConstrS+[DC*(inv(A)-U_s*U_s')-DA*(inv(A)*C*inv(A)-U_s*o_s*U_s') == zeros(size(A))]:'No Spillover';

x_START=[1.0451;

      0.38147;

      0.053452;

      0.15835;

      0.80087;

      75153;

      1.0215e+05;

      68801;

      4964;

      72362;

      -0.2;

      -0.03;

      -0.01;

      -0.6;

      -0.11579;

      -0.0044277;

      -41465;

      -167.12;

      -0.0056313;

      -0.1006;

      -713.07;

      -3746.3];  

 assign([x;y;t],x_START);  

diag_cvx=optimize([s==rvec]+ConstrTRelax+ConstrBounds+ConstrS,s'*s,sdpsettings('solver','gurobi','verbose',1,'usex0',1,'savedebug',1));

[optimizer]

Sorry, wrong copy of last line:

diag_cvx=optimize([s==rvec]+ConstrTRelax+ConstrBounds+ConstrS,s'*s,sdpsettings('solver','fmincon','verbose',1,'usex0',1,'savedebug',1));

Gurobi doesn't support x0...

[Johan Löfberg]

not sure what you mean. among those equallities there are expressions saying x(i)==0, hence lower and upper bounds are 0 (or 5e-16, which means you probably shpuld clean data to ensure you don't have those numbers, as it probably is numerical noise and they should be 0)