Solver not applicable (sedumi)

[woong G. jeong]

Hi.

i've been coding for region of attraction, but i got in trouble  for error :

Warning: Solver not applicable (sedumi)

could u let me know what this problem is ??

Thank you.

% setup model
Ap=-1/3;
Bp=1/3;
Cp=-3;
Kaw=1;
Ki=3;
Kp=5;
A=[Ap-(Bp*Kp*Cp) Bp;-Ki*Cp 0];
B=[-Bp;-Kaw];
F=[-Kp*Cp 1];

n = length(A);
Q=sdpvar(n,n);
v=eye(1,1);
G=sdpvar(1,n);
H=sdpvar(1,n);
gi=G(1,:);
x=sdpvar(n,1);
R=sdpvar(n,n);

opts=sdpsettings('solver','sedumi');
gamma=sdpvar(1);

const1= R*gamma*Q;
M=v*F*Q+(eye(1)-v)*G;
cs1=Q*(A-B*F)'+(A-B*F)*Q+B*M+M'*B';
cs2=[1 gi; gi' Q];
const=const1>=1;
const=[const, Q>0];
const=[const, cs1<0];
const=[const, cs2>=0];
    

% solve
solvesdp(const,gamma,opts);

%Extract numerical solution
gamma=double(gamma);

%plot
plot(x'*double(R)*x<=1,'b');

[Johan Löfberg]

You have a polynomial semidefinite program (R*gamma*Q). Much (much much) harder than linear semidefinite programs which sedumi (and mosek/sdpot3/csdp/dsdp/sdpa...) is built for.

There is basically no good way to solve the problem

You can try

1. A local nonlinear SDP-solver (only option is PENLAB)

2. A global solver (only option is BMIBNB in combination with PENLAB as local solver)

2. Some iterative scheme based on fixing some variables and optimize over others (such as the most simple case of quasi-cconvex problems which can be solved to global optimality by bisection http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Examples.DecayRate). In general, somehow linearize the model and iterate, path-following, trust-regions. Basically, you have to write a nonlinear solver based on repeated linear approximations

[Johan Löfberg]

BTW, strict inequalities has no relevance in YALMIP and associated solvers
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Prepare-your-code

[Johan Löfberg]

and your model is not well-posed. R*Q where R and Q are symmetric matrices, is not necessarily a symmetric matrix, so it is not clear what you mean with R*gamma*Q>=0

...and (eye(1)-v) is zero, so why do you have that term at all?

[woong G. jeong]

Thank you i will try it with non-linear SDP-solver

in fact, constraint is inv(R)<=gamma*Q , (not R*gamma*Q >=1)

it is okay that i write it such as           const1=inv(R)<=gamma*Q

                                                           const=[const1, Q>0]

                                                           const=[const, cs<1]    ?? 

and, i want to solve the problem(not ill-posed problem u mentioned)

how can i make my model well-posed? would u please let me give me any advice?  

[Johan Löfberg]

You should not try to solve this using a generic non linear approach. Instead, you should take the opportunity to learn more about modelling in an SDP setting

To begin with, inv(R)<=gamma*Q should be rewritten using a Schur complement, to arrive at an expression linear in R

For the other constraint, I suspect you can get rid of the nonlinear F*Q terms by performing a standard change of variables where you work with F*Q as a new variable, call it, e.g. Y. It looks like the model will be linear in (Y,Q).

The only remaining problem then would be the product gamma*Q. Unless it is possible to do something clever there, you could minimize gamma by a simple bisection (see decay rate example on wiki)

However, there must be something wrong in your description. R only enters in inv(R)<=gamma*Q as far as I can see, hence trivially it is optimal with R = inf*eye(n) leading to gamma = 0

[woong G. jeong]

Thank you for your advice 

i am trying to change nonlinear to linear using Schur complemet

however, i tried nonlinear using penlab solver (error : Warning: Solver not applicable (penlab))

could u let me know what cause of the error is ? 

[Johan Löfberg]

Your code does not make sense. In the doc-file you sent, you said you wanted the constraint to hold for all v in some set. Now all of a sudden you want to find a binary v such that the constraints hold. Which is it? Moreover, since R only appears in the (1,1) block of a matrix you constrain to be psd, the problem is trivially solved by R = inf*eye(2), gamma = 0+ and Q any matrix satisfying the Lyapunov constraint)

penlab cannot be used with combinatorial (intvar) constraints.

(and don't solve this using penlab. If your intvar construction was a mistake, and you simply want the constraints to hold for v=0 and v=1, then code that, and the only complication is the product gamma*Q, which you can address using bisection (decay-rate example on YALMIP Wiki. But as I said, your model is wrong as it is trivially solved by the solution above)

[dinh nguyenvan]

Hi Johan,

I have a similar problem,

Can I attach the matlab code here and help me to figure out ?

Best Regards,

Dinh

[Johan Löfberg]

Sure

But ask your self if you are multiplying decision variables anywhere. If so, the problem is obviously the fact that you don't have a linear problem

[Sanaz Sabzevari]

Hi,

Thanks for this group,

I wrote the linear problem for optimization, but I do not know wht this code does not worked?

A = [zeros(3) 0.5*eye(3); zeros(3) zeros(3)];

Bv=[2.6912e-5;6.6545e-6;-6.8449e-6];

Bv=Bv/norm(Bv);

Bi=[0 -Bv(3) Bv(2);Bv(3) 0 -Bv(1);-Bv(2) Bv(1) 0];

Biv=[-0.8821e-7;-0.5693e-7;-0.6142e-7];

Biv=Biv/norm(Biv);

Bdoti=[0 -Biv(3) Biv(2);Biv(3) 0 -Biv(1);-Biv(2) Biv(1) 0];

B = [2*Bi eye(3);2*Bdoti 2*Bi];

 

 

Y = sdpvar(6,6);

L = sdpvar(6,6);

 

d=sdpsettings('solver','sdp');

 

% LMI = Y > 0;

% LMI = [LMI, - (Y*A'+Y*L'+A*Y+L*Y') > 0];

LMI = [(Y*A'+Y*L'+A*Y+L*Y') zeros(6,6);zeros(6,6) Y ] ;

Const=[LMI > 0]

solvesdp(Const)

k=double(L)*inv(B);

Could you please help me?

Thanks alot, in advance

[Johan Löfberg]

1. There is no solver called 'sdp'

https://yalmip.github.io/allsolvers/

2. There is no command today called solvesdp

https://yalmip.github.io/tutorial/basics/

3. Strict inequalities are not supported

4. You don't want L to be symmetric (you've defined it as symmteric now)

5. Weird to call ithe expression "LMI" when you obviously have a BMI (which you cannot solve) as you multiply Y and L 

6. You're defining matrix B which you never use

[Johan Löfberg]

7. To define X= 0 and Y>=0 you simply write [X>=0,Y>=0], Don't create a large matrix to hold both variables

[Sanaz Sabzevari]

Thanks a lot for your response, 

I made some changes but k is zero

clear all

A = [zeros(3) 0.5*eye(3); zeros(3) zeros(3)];

Bv=[2.6912e-5;6.6545e-6;-6.8449e-6];

Bv=Bv/norm(Bv);

Bi=[0 -Bv(3) Bv(2);Bv(3) 0 -Bv(1);-Bv(2) Bv(1) 0];

Biv=[-0.8821e-7;-0.5693e-7;-0.6142e-7];

Biv=Biv/norm(Biv);

Bdoti=[0 -Biv(3) Biv(2);Biv(3) 0 -Biv(1);-Biv(2) Biv(1) 0];

B = [2*Bi eye(3);2*Bdoti 2*Bi];

 

 

Y = sdpvar(6,6,'full');

L = B*Y;

k=sdpvar(6,6,'diagonal','real');

 

ops = sdpsettings('solver','sedumi');

 

% LMI = Y > 0;

% LMI = [LMI, - (Y*A'+L'*k'+A*Y+k*L) > 0];

LMI = [- (Y*A'+L'*k'+A*Y+k*L) zeros(6,6);zeros(6,6) -Y ] ;

Const=[LMI <0]

solvesdp(Const,3)

value(k)

value(Y)

value(L)

[Johan Löfberg]

Looks very wrong. Isn't this a standard stabilization problem, i.. you want to find P>0 and K such that P(A+BK) + (A+BK)P < 0

Congruence with Y = P^-1 leads to AY + BKY + (AY+BKY)' < 0

Define KY = L leading to the (AY + BL) + (AY+BL)' < 0 which is an LMI in the decision variables Y and L. Accept that you cannot have strict inequalities and replace with (AY + BL) + (AY+BL)' <= 0., Y>=0. Understand that this problem is homogeneous and thus Y=L=0 is a solution, but that can easily be fixed by simply using Y>=eye(n) instead

You are currently defining what I see as a Lyapunov function Y to be non-symmetric. That makes no sense. 

Saying k is 'real' is redundant, as 'real' is default.  But generally I think you're doing something wrong in defining a diagonal matrix k and defining your constraint, as you still have a nonconvex nonlinear constraint as you have products between k and Y. You still don't seem to appreciate what the L in LMI stands for...

some more weirdness

1. You minimize the objective "3" ??

2. You define an options structure but you never use it

3. You're still using solvesdp

4. You're still using strict inequalities

5. You're still unnecessarily blocking variables to define completely separate constraints

Did you really read the tutorials?

[Sanaz Sabzevari]

Ohhh great Thanks for your comments,

Yes, Right now it works and I do all the comments :)

Many Thanks for creating this group and your fast and efficient response

[Sanaz Sabzevari]

Just, sorry could you help me how could I use moseq as a solver?

Because I install moseq and I do not know why I have this error :

Unkown problems in solver (turn on 'debug' in sdpsettings)

Thanks a bunch :)

[Sanaz Sabzevari]

addpath(genpath('C:\Program Files\Mosek'))

clear all

A = [zeros(3) 0.5*eye(3); zeros(3) zeros(3)];

%Bv=[2.6912e-5;6.6545e-6;-6.8449e-6];

Bv=[2.535862e-5;6.9330e-6;2.1018e-5];

Bv=Bv/norm(Bv);

Bi=[0 -Bv(3) Bv(2);Bv(3) 0 -Bv(1);-Bv(2) Bv(1) 0];

%Biv=[-0.8821e-7;-0.5693e-7;-0.6142e-7];

Biv=[-1.0426e-7;0.01304e-7;5.6791e-7];

Biv=Biv/norm(Biv);

Bdoti=[0 -Biv(3) Biv(2);Biv(3) 0 -Biv(1);-Biv(2) Biv(1) 0];

C = [2*Bi eye(3);2*Bdoti 2*Bi];

 

 

P = sdpvar(6,6,'full');

L = sdpvar(6,6,'full');

 

d=sdpsettings('solver','mosek');

 

LMI = P -1e-6 >= 0;

LMI = [LMI,  (A'*P+C'*L'+P*A+L*C)-1e-6  <= 0];

%LMI = [(Y*A'+Y*L'+A*Y+L*Y') zeros(6,6);zeros(6,6) Y ] ;

%Const=[LMI > 0]

solvesdp(LMI)

k=value(L)*inv(value(P));

eig(double(P))

eig(double(k))

This is the new code that I wrote.

Super thanks, in advance.

[Johan Löfberg]

So have you confirmed that mosek actually work on its own? Run the command mosekdiag and mosek will check its installation

[Johan Löfberg]

1. You're still using solvesdp. Look at the eamples, no such command is used. It is calle optimize

2. P >= 1e-6 is equivalent to P >= 1e-6*ones(n) which is not what you want. You want  P >= 1e-6*eye(n)

[Johan Löfberg]

...and since P is a Lyapunov matrix, it should not be a full matrix.

[Sanaz Sabzevari]

I applied all of your comments, but still it does not worked 

clear all

A = [zeros(3) 0.5*eye(3); zeros(3) zeros(3)];

%Bv=[2.6912e-5;6.6545e-6;-6.8449e-6];

Bv=[2.535862e-5;6.9330e-6;2.1018e-5];

Bv=Bv/norm(Bv);

Bi=[0 -Bv(3) Bv(2);Bv(3) 0 -Bv(1);-Bv(2) Bv(1) 0];

%Biv=[-0.8821e-7;-0.5693e-7;-0.6142e-7];

Biv=[-1.0426e-7;0.01304e-7;5.6791e-7];

Biv=Biv/norm(Biv);

Bdoti=[0 -Biv(3) Biv(2);Biv(3) 0 -Biv(1);-Biv(2) Biv(1) 0];

C = [2*Bi eye(3);2*Bdoti 2*Bi];

 

 

P = sdpvar(6,6);

L = sdpvar(6,6,'full');

 

d=sdpsettings('solver','sedumi');

 

LMI = P >= 1e-8*eye(6) ;

LMI = [LMI, -(A'*P+C'*L'+P*A+L*C)  >= 1e-8*eye(6)];

%LMI = [(Y*A'+Y*L'+A*Y+L*Y') zeros(6,6);zeros(6,6) Y ] ;

%Const=[LMI > 0]

optimize(LMI)

k=value(L)*inv(value(P));

eig(double(P))

eig(double(k))

This is a new code

Super thanks for your assistant 

[Sanaz Sabzevari]

When I used this command 'mosekdiag' this message pops up:

Matlab version: 8.5.0.197613 (R2015a)

Architecture  : PCWIN64

The mosek optimizer executed successfully from the command line:

 

MOSEK Version 8.1.0.42 (Build date: 2018-3-13 02:09:02)

Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com

Platform: Windows/64-X86

 

FlexLM

Version                  : 11.14

Hostname                 : AAU109226

Host ID                  : "0050b6ce319a 4c3488318c4d 4c3488318c4c 4c3488318c50"

License path             : C:\Users\sanaz\mosek\mosek.lic

 

Operating system variables

PATH                     :

 

*** No input file specfied. No optimization is performed.

 

Return code - 0  [MSK_RES_OK]

 

mosekopt: C:\Program Files\Mosek\8\toolbox\r2014a\mosekopt.mexw64

Error using mosekdiag (line 58)

mosekopt returned an error: 1008. [MSK_RES_ERR_MISSING_LICENSE_FILE]

[Sanaz Sabzevari]

and of course with moseq solver there is still this problem

[Johan Löfberg]

Hence you know what the problem is

[Johan Löfberg]

works here

you're still not using the options structure you create though

[Sanaz Sabzevari]

I do not understand what you mean about this?

[Sanaz Sabzevari]

yalmiptime: NaN

    solvertime: NaN

          info: 'Unknown problem in solver (Turn on 'debug' in sdpsettings) (Reference to non-existent field 'sol'.)'

       problem: 9

This message pops up here after running the problem with both sedumi and mosek solvers 

[Johan Löfberg]

You're not using the options, you're only creating it.

https://yalmip.github.io/tutorial/basics/

[Johan Löfberg]

1 You've already run mosekdiag and been told what the problem is

2. You only create an options structure defining the solver choice, but you never use it in optimize, hence YALMIP uses the best solver it can find, which is (your failed installation of) mosek

[Johan Löfberg]

close, but you did not read the tutorial or manual carefully

https://yalmip.github.io/tutorial/basics/

https://yalmip.github.io/command/optimize/

and it's 'mosek'

[Johan Löfberg]

mosek support sounds like a reasonable start

[Sanaz Sabzevari]

yes, but with Sedumi I have a result and for mosek it is infinitive

[Sanaz Sabzevari]

I do not know what should I do then?

Why this is a very different response and may I trust to the sedumi solver?

[Johan Löfberg]

I have no idea what you're saying.

If you want to use mosek, you contact mosek support to find out why the license isn't found.

Very different response?

[Sanaz Sabzevari]

yes sedumi has a response after 4th iteration but mosek it ia NAN.

Just, I do not know that the license has a prblem or other thing ...

[Johan Löfberg]

you haven't solved any problem with mosek, it detects a license problem

[Sanaz Sabzevari]

But I install trial license for 30 days and this is a 3rd day.

[Johan Löfberg]

As I said, this is an issue for mosek support.

[Sanaz Sabzevari]

Super thanks for all of your assistance. 

[Sanaz Sabzevari]

Hi,

Is it possible to have a scheduling parameter p variable with respect to time and apply it in Yalmip? Thereby, we have the solution dependant on time?

I mean could we get the solution functioning with time in Yalmip?

Best

[Johan Löfberg]

No.

You have to derive the optimization problem, i.e. derive suitable theory to derive a finite-dimensional problem, similar to this

https://yalmip.github.io/example/lpvstatefeedback/

[Sanaz Sabzevari]

Thanks a lot.

[Sanaz Sabzevari]

Hi again,

For LMI that the bound is >0, we must use >=1, and when we change this bound to >=1e-3 the response is very different. So, what criteria is important to choose this?

Super thanks 

[Johan Löfberg]

Your question is " How can I replace a strict inequality > 0 with a non-strict >=". 

There is no simple way. If the problem is homogeneous, any margin works as it cannot influence feasibility, hence >= 1  (>= eye(n) for semidefnite constraints) is the obvious choice. For non-homogeneous problems, the margin will affect the feasible region and thus it has to be selected with care. Too large and it will influence the feasible region and might render the problem infeasible or result in sub-optimal solution, but too small and it will drown in the general tolerances of the solvers and you effectively have a non-strict condition

[rp1]

Hi Johan,

I am having a similar problem in trying to find an upper bound for a series of differential equations of a large order system. For now I have a set differential system with 9 variables which uses two functions to reduce the system to the first 4 variables (4 is an input but this number can change) plus a residual. The result is a new system of four variables a_1,...a_4 and the residual q^2 for which I am trying to find a bound on the energy of the system equal to the sum of (a_n)^2 + q^2.

 I am also getting an error saying that the solver is non-applicable, however the introduction of a z variable which is done towards the end of the script is such that to make the problem linear in the decision variables a and qsquared as given below. If you could give any insight into what the issue may be that would be great help. I have attached my script below, many thanks.

clear, clc, clear yalmip

 

opts = sdpsettings();

opts.solver = 'sdtp';

opts.cdcs.maxIter = 1e4;

opts.cdcs.solver = 'sos';

 

Re = 400; % CHANGE Re TO DESIRED VALUE

 

% Generating the full system

[N_ijk, L_ij, B_i] = Shear_Model_9_Mode(Re);

Size_N_ijk = size(N_ijk);

N = Size_N_ijk(1);

 

% Use UODESys tool to undergo variable transformation

Variable_Transformation

 

% Use UODESys tool to find bound for Theta and value for KAPPA

Bounds_Calculation

 

% Define column vector of variables [a_1; ... ;a_n]

a = sdpvar(n,1);

 

% Represent the 3D array N_hat_xyz as a cell array. 

N_hat_xyz_final_array = cell(1,n);

for x = 1:n

    N_hat_xyz_final_array{1,x} = zeros(n,n);

    for y = 1:n

        N_hat_xyz_final_array{1,x}(y,:) = N_hat_xyz_final(x,y,:);

    end

end

 

% Represent N_hat_xyz_final_array as a 2D matrix.

nnl_mat = [N_hat_xyz_final_array{1}(:)'; N_hat_xyz_final_array{2}(:)'; N_hat_xyz_final_array{3}(:)'; N_hat_xyz_final_array{4}(:)'];

 

% Define f=da/dt for new truncated system using Transformed N,L and B

% matrices

 

a2 = a*a';

lin = B_hat_x_final + L_hat_xy_final*a;

nnl = nnl_mat*a2(:);

f = lin + nnl;

 

% Assign values to the sdpvar decision variables.

assign(c_1,value(c_1))

assign(c_2,value(c_2))

assign(c_3,value(c_3))

 

 

sdpvar r;

vector_of_variables = [a;r];

 

%Create storage function V that is a function of u and q^2

[V_6,c_6,mono_6] = polynomial([vector_of_variables],6,1);

 

%Jacobian of V

F = jacobian(V_6,[vector_of_variables]); 

 

 

%Assign q^2 as variable, equal to norm(v)/2

sdpvar q;

qsquared = q^2;

 

%Define phi as sum of squares of u1,...un + q^2

phi = norm_u_sq + qsquared;

 

%replace r with qsquared in jacobian F

replace(F,r,qsquared);

 

%Assign variable z0 and z

sdpvar z0

z = sdpvar(n,1);

 

%Assing expression for P

norm_a_sq = a'*a;

P = c_1*qsquared + c_2*norm_a_sq*2*qsquared + c_3*qsquared^2*4;

 

 

 

 

%Solve for upperbound

sdpvar upperbound

Constupper=sos(-(((P*z0^2)+z'*z)*(F(1:n)*f + (F(n+1)*KAPPA*qsquared) + phi - upperbound)') - (2*P*z0*(F(1:n)-F(n+1)*a')*z));

solvesos(Constupper, upperbound,[opts],[a;qsquared;upperbound;c_6]);

[Johan Löfberg]

Your code makes no sense. There is no solver called 'sdtp', and even if it was, why are you then editing options in the solver 'cdcs', 

Since you haven't supplied reproducible code, it is impossible to investigate further, although at a first glance the code looks weird (you assign values to c_1 before you solve any problem so value(c_1) will be nan, and the assigned values will not be used anyway since you haven't told the solver to use initial guesses, and the solver you want to use (cdcs?)  does not support initial guesses anyway...)

[rp1]

Hi Apologies, that should say 'cdcs' instead of sdtp on that line, as the system is using the solver cdcs.

I should have explained that the beginning two functions: Variable_transformation and bounds_calculation produce the arrays N_hat_xyz_final_array, L_hat_xy_final_array and B_hat_x_final_array, as well as the values for c_1,..c_3 so it can be assumed thats assign(c_1,value(c_1)) generates the output from the function bounds_calculation which gives c_1 a constant value, and the same goes for c_2, c_3.

[Johan Löfberg]

Why are you even using the commands assign. Those lines are are redundant in your code

[rp1]

Ahh you're right I can remove them thank you. 

Do you have any clue as to why why I get the error solver not applicable (cdcs) however, debugging the script everything runs fine untill the final line of the code however I am unsure as to why the equation set up may be invalid.

[Johan Löfberg]

one thing is that your polynomial is nonlinear in the parameter qsquared, hence the resulting SDP will be nonlinear. Same thing with a.

[rp1]

Okay thank you for your help. Much appreciated

[Johan Löfberg]

it looks like you are confused with the idea of parameters in the sos problem, i.e. that you have confused about what are parameters (to be found) and independents (which the sos decomposition should be performed over) in the sos problem