Solve FMO - a nonlinear SDP problem by yalmip and penlab

[jingqiao hu]

Dear professor:

I have two questions about nonlinear SDP optimization.

1. I'm solving FMO(free material optimization) problem. It's a nonlinear SDP problem and the formulation is:

I used the PENLAB to solve. However, the iteration output is 

************* Start of outer step   1 **********

object(x_  0) =                      Inf

||grad(x)||_2 =   4.2788351142571500E+01

     ----

Chol fact: OK, pert=1.0000e-06

LS (pen): -1.0000e+06, 20 steps, Step width: 35.200000

object(x_  1) =                      Inf

||grad(x)||_2 =   1.8333380466026046E+01

     ----

Chol fact failure (312), last pert: 5.629500e+08, giving up

FAILURE: Newton solver (0) cannot continue (flag 1)

*****!!!!!!! Result of outer step   1 !!!!!*****

Objective                -3.5200000000000000E+07

Augmented Lagrangian                         Inf

|f(x) - f(x_old)|         3.5200000000000000E+07

|f(x) - Lagr(x)|                             Inf

Grad augm. lagr.          1.8333380466026046E+01

Feasibility (max)         1.0000000512000000E+00

Feasibility eqx      

Feasibility ineq          1.0000000000000000E-04

Feasibility box      

Feasibility m.ineq        1.0000000512000000E+00

Complementarity           2.5000000000000000E-01

Minimal penalty           7.6007250103963326E-52

Newton steps                                   1

Inner steps                                   53

Linesearch steps                              20

Time of the minimization step          36.7969 s

  - factorizations in the step         3.40625 s

*****!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*****

 

Warning: Unconstr/eqconstr. min failed (flag 2), remaing attepmts: 2

Is it possible that the problem is too large or the derivatives that yalmip deduced is wrong?

2. The corresponding of linear formulation of FMO problem is 

I know it has very large-scale matrix inequalities which are difficult to deal with in computations.  But I still try to optimize it using sdpt3 and the result is different compared with answer.

Is it possible that the matrix inequalities is too difficult to deal with so that something wrong in the computation?

Thanks in advance.

Best regards.

[Johan Löfberg]

I've never got penlab to work well, it is very unstable

Regarding wrong answer, you would have to define better what wrong means

[jingqiao hu]

Thanks for reply. The details of wrong answer are as follows:

I drew the colormap of the trace of variables E. Area in red circle behaves inappropriate. There should be some material(yellow area) in the red circle.

The following figure is the right answer. The corresponding color in the two colormap is different but the whole color gamut is same.

And what about pennon in tomlab? Is it stable?

[jingqiao hu]

Could you recommend some solver to solve nonlinear SDP problem? Thank you.

[Johan Löfberg]

Well, the question is more what is wrong with the solution in the context of the optimization process, i.e. what does check return when you check the constraints, what does the solver display, and what diagnostics code are returned. If all those look ok, it is a sign that your model is wrong

[Johan Löfberg]

No, nonlinear semidefinite is extremely hard and there is no generic solver which works as stable as you would have for linear SDP. Kocvara has made a lot of work on this for decades, but it still doesn't solve everything out of the box (as you've noticed), and it is still a research problem

[jingqiao hu]

Thanks for the explanation. The linear SDP formulation in the question 2 is solved by SDPT3. And the information presented is 

*******************************************************************

   SDPT3: Infeasible path-following algorithms

*******************************************************************

 version  predcorr  gam  expon  scale_data

   HKM      1      0.000   1        0    

Is it possible that the Infeasible path-following algorithms bring numerical instability in sdpt3?

[Johan Löfberg]

If that's all you see, the solver has crashed, and you should have an error code in the diagnostics.

Try mosek, it is better typically

[jingqiao hu]

The attachment is whole output. The solver isn't crashed and it finished the optimization. 

I have double check my model and I think it's right. Thus I think the problem is the solver itself.

Thanks for reply and I will try mosek.

[Johan Löfberg]

sdpt3 has terminated without convergence

As I said, try mosek instead

[jingqiao hu]

The results of formulation in the question 2 solved by MOSEK and SDPT3 both are DUAL INFEASIBILITY.

The following code is the corresponding matlab code. 

The matrix in constraint [gamma, f^T; f, A(E)]>=0 is huge. Is it possible that semidefinite constraint for large matrix lead to dual infeasibility?

%%********************************************************************

% linear SDP but slow and huge

% min(E) sum(trace(E_i))

% s.t.  Ei>=0

%       l<= trace(E_i) <=u

%       [gamma, f^T; f, A(E)]>=0

%%********************************************************************

dbstop if error

clear all

yalmip('clear') 

nelx = 20;

nely = 10;

m = nelx*nely;

l = 1e-4;

u = 1;

gamma = 170;

% ======================= FEA ==========================================

fixeddofs = 1:2*(nely+1);

alldofs = 2*(nely+1)*(nelx+1);

freedofs = setdiff(1:alldofs,fixeddofs);

% loadnid = 2*(nely+1)*(nelx+1);

loadnid = (nely+1)*nelx+1+(nely+1)*(nelx+1);

F = sparse(loadnid,1,-1,alldofs,1);

% for assmble 

nodenrs = reshape(1:(1+nelx)*(1+nely),1+nely,1+nelx);

edofVec = reshape(2*nodenrs(1:end-1,1:end-1)+1,m,1);

edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],m,1);

iK = reshape(kron(edofMat,ones(8,1)).',64*m,1);

jK = reshape(kron(edofMat,ones(1,8)).',64*m,1);

%============================== FMO ===============================

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

ops.showprogress = 1; %1为设置显示yalmip现在在做什么

ops.verbose = 2; %设置显示信息程度，1为适度显示，2为完全显示。

ops.debug = 1;

trE = 0;

cons1 = [];

cons2 = [];

for i = 1:m    

    E{i} = sdpvar(3); % 3*3 symmetric

    

    assign(E{i},repmat(0.1*u,3,3));

    cons2 = [cons2,E{i} >= 0]; % constraint 2: Ei>=0        

    cons1 = [cons1,l<= trace(E{i}) <=u]; % constraint 1: l<= trace(E_i) <=u

    

    trE = trE+trace(E{i});

    

    tempE = E{i};

    sE(:,i) = [tempE(1,1);tempE(1,2);tempE(1,3);tempE(2,2);tempE(2,3);tempE(3,3)];

end

sK = getK(sE);

K = sparse(iK,jK,sK);

K = (K+K')/2;

% constraint 3: [gamma, f^T; f, A(E)]>=0

cons3 = [gamma,F';F,K] >= 0; 

obj = trE;

sol = optimize([cons1,cons2,cons3],obj,ops);

[Johan Löfberg]

Doesn't look like a particularly large problem according to the log

Undefined function or variable 'getK'.

 

[jingqiao hu]

Sorry for the omission. Here is getK() function.

[Johan Löfberg]

It's simply not feasible

Let gamma be a parameter, and simply solve feasibility problem

>> sdpvar gamma

>> optimize([cons1,cons2,[gamma F';F K]>=0])

Note that there are a bunch of negative eigenvalues, i.e. despite simply solving a feasibility problem, which typically leads to interior solutions, it has ended up with a slightly infeasible

>> eig(value([gamma F';F K]))

ans =

    -8.965346389965378e-09

    -2.197306043384670e-12

    -1.411122297642595e-13

     2.824596349326144e-03

     5.866168694096923e-03
...

my guess is that there is a structural problem with the model causing it to have some singular features

try to minimize gamma to see how far of the value 170 is

optimize([cons1,cons2,[gamma F';F K]>=0],gamma)

and it runs into numerical problems, and terminates with a solution with gamma = 7179, and also here 3 problematic eigenvalues 

>> eig(value([gamma F';F K]))

ans =

    -1.895104879418130e-07

    -3.988447448666848e-14

     7.432205916707894e-13

     2.355337978541212e-03

Despite setting gamma to infinity essentially, you will still not satisfy the original constraint, so there appears to be something structurally wrong

>> assign(gamma,1e32);

>> min(eig(value([gamma F';F K])))

ans =

    -3.053257919279714e+16

[jingqiao hu]

Thank you so much for your detailed explanation, Johan. It helps me a lot. I will double check my model.