Nonlinear problem using YALMIP+IPOPT

Guori Huang

unread,

Apr 1, 2016, 7:40:14 AM4/1/16

to YALMIP

When I use YALMIP+IPOPT to solve a large-scale nonlinear problem, it takes too much time to solve it, exceed CPU_max_time. How can I improve computing efficiency?

Johan Löfberg

unread,

Apr 1, 2016, 7:41:23 AM4/1/16

to YALMIP

Impossible to answer. Create a smaller less complicated model...

Guori Huang

unread,

Apr 2, 2016, 3:28:20 AM4/2/16

to YALMIP

I am a new learner, when I added the following constraints:

F_op = F_op + [ S_pipe_op == sign( sparse(A_connect_gasnode_pipe_TR) * pressure) .* sqrt( S_pipe_op_square ) ];

F_op = F_op + [ S_pipe_op_square == k_pipe.^2 .* sign( sparse(A_connect_gasnode_pipe_TR) * pressure) .*...

( sparse(A_connect_gasnode_pipe_TR ) * pressure )...

.* ( sparse( A_connect_gasnode_pipe_TR_sum ) * pressure) ];

in that ,

A_connect_gasnode_pipe_TR is a 22*17 constant matrix;

pressure = sdpvar(17,10*3);

S_pipe_op_square = sdpvar(22,10*3);

S_pipe_op = sdpvar(22,10*3);

k_pipe is a 22*(10*3) constant matrix;

A_connect_gasnode_pipe_TR_sum is a 22*17 constant matrix;

warning:

ans =

Model creation failed (Model not available in constraint #1033 at level 1)

ans =

14

ans =

Model creation failed

what is meannings of the warning?

Thank you very much.

Johan Löfberg

unread,

Apr 2, 2016, 3:40:30 AM4/2/16

to YALMIP

sqrt is convexity aware, and since you are using it in a nonconvex setting, it fails

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Extra.SQRT

Guori Huang

unread,

Apr 2, 2016, 7:09:56 AM4/2/16

to YALMIP

Thank you very much! Your advice is very helpful for me. If I want to get dual multipliers of the above constraints, how should I do? because I need to return the dual multipliers to major program. Do you have any advice? Thank you!

Johan Löfberg

unread,

Apr 2, 2016, 10:02:16 AM4/2/16

to YALMIP

the command dual

Guori Huang

unread,

Apr 3, 2016, 2:05:21 AM4/3/16

to YALMIP

Is your meaning that dual variable is equal to dual multiplier? both of them is equivalent.

Johan Löfberg

unread,

Apr 4, 2016, 1:55:32 AM4/4/16

to YALMIP

dual variable = dual multiplier = lagrange multiplier etc. Different names

Guori Huang

unread,

Apr 4, 2016, 9:07:03 AM4/4/16

to YALMIP

%% S_DNGFP_op is a 6*30 variable matrix; S_DNGFP_max is a 6*30 constant matrix

F_op = F_op + [( S_DNGFP_op <= S_DNGFP_max ) : 'DNGFP'];

%% S_GasStation_op is a 2*30 variable matrix; S_GasStation_max is a 2*30 constant matrix

F_op = F_op + [( S_GasStation_op <= S_GasStation_max) : 'GasStation'];

%% A_connect_node_DNGFP is a 54*6 constant matrix; P_DNGFP_op and Q_DNGFP_op are 6*30 variable matrixs;

%% A_connect_node_grid is a 54*2 constant matrix; P_grid_op and Q_grid_op are 2*30 variable matrixs;

%% A_connect_node_line is a 54*61 constant matrix; P_line_op and Q_line_op are 61*30 variable matrixs;

%% P_d and Q_d are 54*30 constant matrixs; P_d_cut_op, P_d_analog1 and P_d_analog2 are 54*30 variable matrixs;

%% Q_d_cut_op, Q_d_analog1 and Q_d_analog2 are 54*30 variable matrixs;

F_op_sub_analog = [(sparse( A_connect_node_DNGFP ) * P_DNGFP_op...

+ sparse( A_connect_node_grid ) * P_grid_op...

- sparse( A_connect_node_line ) * P_line_op...

- P_d + P_d_cut_op + P_d_analog1 -P_d_analog2 == 0 ): 'elePflow' ];

F_op_sub_analog = F_op_sub_analog + [ ( sparse( A_connect_node_DNGFP ) * Q_DNGFP_op...

+ sparse( A_connect_node_grid ) * Q_grid_op...

- sparse( A_connect_node_line ) * Q_line_op - Q_d...

+ Q_d_cut_op + Q_d_analog1 -Q_d_analog2 == 0 ): 'eleQflow' ];

%% A_connect_gasnode_GasStation is a 17*2 constant matrix; A_connect_gasnode_DNGFP is a 17*6 constant matrix;

%% A_connect_gasnode_pipe is a 17*22 constant matrix; S_d_g is a 17*30 constant matrix;

%% S_pipe_op is a 22*30 variable matrix; S_d_g_analog1 and S_d_g_analog2 are 17*30 variable matrixs;

F_op_sub_analog = F_op_sub_analog + [ (sparse(A_connect_gasnode_GasStation) * S_GasStation_op...

- sparse(A_connect_gasnode_DNGFP) * S_DNGFP_op/y_efficiency...

- S_d_g - sparse(A_connect_gasnode_pipe) * S_pipe_op...

+ S_d_g_analog1 - S_d_g_analog2 ==0 ): 'gasflow' ];

DNGFP = dual(F_op('DNGFP'));

GasStation = dual(F_op('GasStation'));

elePflow = dual(F_op_analog('elePflow'));

eleQflow = dual(F_op_analog('eleQflow'));

gasflow = dual(F_op_analog('gasflow'));

Results of DNGFP,GasStation,elePflow,eleQflow and gasflow are NaN. What is wrong with that?

Johan Löfberg

unread,

Apr 4, 2016, 9:16:38 AM4/4/16

to YALMIP

Misinterpreted you. YALMIP does not recover duals from nonlinear programs solved using ipopt.

Guori Huang

unread,

Apr 4, 2016, 9:40:05 AM4/4/16

to YALMIP

Is there any other way to get duals？

Johan Löfberg

unread,

Apr 4, 2016, 9:44:13 AM4/4/16

to YALMIP

You can always get the numerical values by having YALMIP catching all data returned from the solver by using the option savesolveroutput, and then look at the data in the struct returned by optimize

Johan Löfberg

unread,

Apr 4, 2016, 9:53:55 AM4/4/16

to YALMIP

BTW, it is completely redundant to do sparse(A)*x in any YALMIP code, as everything is saved as sparse internally in YALMIP anyway.

Guori Huang

unread,

Apr 4, 2016, 9:48:10 PM4/4/16

to YALMIP

Is there any other solvers which is able to recover duals from nonlinear programs?

Guori Huang

unread,

Apr 5, 2016, 12:58:10 AM4/5/16

to YALMIP

And how to use the option savesolveroutput?

Guori Huang

unread,

Apr 8, 2016, 3:22:58 AM4/8/16

to YALMIP

Could you please help me to solve the problem?

Johan Löfberg

unread,

Apr 8, 2016, 4:09:36 AM4/8/16

to YALMIP

Help with what how? As I said, YALMIP does not recover duals for nonlinear programs. You can extract the numerical values directly from the solver using the savesolveroutput, but that's all you have.

Guori Huang

unread,

Apr 8, 2016, 4:30:00 AM4/8/16

to YALMIP

When I use the command 'sol_op_analog.solveroutput.info.lambda', answers are zero.

Guori Huang

unread,

Apr 8, 2016, 4:30:48 AM4/8/16

to YALMIP

Is that the correct duals?

Johan Löfberg

unread,

Apr 8, 2016, 4:32:21 AM4/8/16

to YALMIP

Those are the values that ipopt returned. Whether they are correct I don't know of course.

Guori Huang

unread,

Apr 8, 2016, 4:33:56 AM4/8/16

to YALMIP

Anyway, thank you very much.

Guori Huang

unread,

Apr 8, 2016, 6:51:59 AM4/8/16

to YALMIP

By solveroutput, I find the duals, but how to determine which constraints the multiplier corresponds to？

Johan Löfberg

unread,

Apr 8, 2016, 7:03:52 AM4/8/16

to YALMIP

That's the problem. YALMIP does not recover them for you so you have to hope that the duals are in the same order as your constraints, and that YALMIP hasn't introduced any auxilliary modelling constraints etc. Short story: No support for duals in nonlinear programs.

Guori Huang

unread,

Apr 8, 2016, 9:07:54 AM4/8/16

to YALMIP

According to the situation, I prepare for linearizing nonlinear constraints, however I couldn't linearizing the square term of the objective function. In this situation, are there suitable solvers to get the duals?

Johan Löfberg

unread,

Apr 9, 2016, 2:47:42 AM4/9/16

to YALMIP

Couldn't linearize?

Anyway, YALMIP recovers duals from all QP solvers

Message has been deleted

Guori Huang

unread,

Apr 11, 2016, 4:00:07 AM4/11/16

to YALMIP

Yes.

Whether Yalmip can execute a variable judgment statement?

for example,

if a<b

F=F+[ t == a-b ];

end

where a,b,t are variables

Johan Löfberg

unread,

Apr 11, 2016, 4:18:40 AM4/11/16

to YALMIP

You will ahve to introduce binary variables and do logic programming via big-M models. See the wiki

E.g implies(a<=b, t == a-b), which manually means you have a binary d and [implies(a<=b,d), implies(d, t == a-b)], which then can be translated to big-M models

Guori Huang

unread,

Apr 11, 2016, 4:19:46 AM4/11/16

to YALMIP

Could YALMIP be used in this way to achieve it?

a = sdpvar(5,5);

b = sdpvar(5,5);

tt= sdpvar(5,5);

F= [ iff(a>=b, tt==a-b) ];

Johan Löfberg

unread,

Apr 11, 2016, 4:25:32 AM4/11/16

to YALMIP

Yes, but remembering that these things are sensitive to tolerances etc, and requires tight explicit bounds on all involved variables for the big-M model to be somewhat tight.

Guori Huang

unread,

Apr 12, 2016, 6:59:28 AM4/12/16

to YALMIP

Is there any way to check error by YALMIP?

E.g.

Tried aggregator 2 times.

MIP Presolve eliminated 4165 rows and 2554 columns.

MIP Presolve modified 1268 coefficients.

Aggregator did 1291 substitutions.

Reduced MIP has 1886 rows, 1236 columns, and 15005 nonzeros.

Reduced MIP has 600 binaries, 0 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.02 sec. (16.08 ticks)

Probing fixed 0 vars, tightened 28 bounds.

Probing time = 0.00 sec. (1.49 ticks)

Tried aggregator 2 times.

MIP Presolve eliminated 8 rows and 1 columns.

MIP Presolve modified 415 coefficients.

Aggregator did 1 substitutions.

Reduced MIP has 1877 rows, 1234 columns, and 14983 nonzeros.

Reduced MIP has 598 binaries, 0 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.01 sec. (7.09 ticks)

Probing time = 0.00 sec. (0.62 ticks)

Clique table members: 191.

MIP emphasis: balance optimality and feasibility.

MIP search method: dynamic search.

Parallel mode: deterministic, using up to 4 threads.

Root relaxation solution time = 0.02 sec. (7.39 ticks)

Nodes Cuts/

Node Left Objective IInf Best Integer Best Bound ItCnt Gap

0 0 infeasible 367

Root node processing (before b&c):

Real time = 0.05 sec. (11.88 ticks)

Parallel b&c, 4 threads:

Real time = 0.00 sec. (0.00 ticks)

Sync time (average) = 0.00 sec.

Wait time (average) = 0.00 sec.

------------

Total (root+branch&cut) = 0.05 sec. (11.88 ticks)

Tried aggregator 3 times.

MIP Presolve eliminated 5293 rows and 2915 columns.

MIP Presolve modified 1314 coefficients.

Aggregator did 1423 substitutions.

Reduced MIP has 626 rows, 743 columns, and 5111 nonzeros.

Reduced MIP has 120 binaries, 0 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.02 sec. (17.90 ticks)

Probing fixed 0 vars, tightened 16 bounds.

Probing time = 0.02 sec. (0.22 ticks)

Tried aggregator 1 time.

Presolve time = 0.00 sec. (1.76 ticks)

Probing time = 0.00 sec. (0.22 ticks)

Clique table members: 60.

MIP emphasis: balance optimality and feasibility.

MIP search method: dynamic search.

Parallel mode: deterministic, using up to 4 threads.

Root relaxation solution time = 0.00 sec. (2.04 ticks)

Nodes Cuts/

Node Left Objective IInf Best Integer Best Bound ItCnt Gap

0 0 infeasible 142

Root node processing (before b&c):

Real time = 0.02 sec. (3.31 ticks)

Parallel b&c, 4 threads:

Real time = 0.00 sec. (0.00 ticks)

Sync time (average) = 0.00 sec.

Wait time (average) = 0.00 sec.

------------

Total (root+branch&cut) = 0.02 sec. (3.31 ticks)

Hmm,something went wrong!

ans =

Infeasible problem (CPLEX-IBM)

ans =

1

There is no useful error messages.

Johan Löfberg

unread,

Apr 12, 2016, 7:34:11 AM4/12/16

to YALMIP

What do you mean. There is a very clear error mesage telling you that the model is infeasible.

Guori Huang

unread,

Apr 12, 2016, 7:44:39 AM4/12/16

to YALMIP

I mean what are reasons resulting that the model is infeasible.

Johan Löfberg

unread,

Apr 12, 2016, 8:27:26 AM4/12/16

to YALMIP

The standard way to debug this is to start removing constraints from the model, until you find a small amount of constraints which together are infeasible. From there, think hard and try to figure out where you've gone wrong with those. Alternatively, add slacks on constraints and minimize hte sum of these slacks. If only a few of the slacks are non-zero, you most likely found the offending constraint for further analysis

Guori Huang

unread,

Apr 13, 2016, 2:33:17 AM4/13/16

to YALMIP

Warning: You have unbounded variables in an implication leading to a lousy big-M relaxation.

But I think all involved variables have explicit bound constraints.

Procedures are as follows:

%% +++++变量

S_GasStation_op = sdpvar(2,10*3); %门站的出力

P_DNGFP_op = sdpvar(6,10*3); %DNGFP的出力P

P_grid_op = sdpvar(2,10*3);% 外部电网出力P

P_d_cut_op = sdpvar(54,10*3); %节点切负荷P

P_line_op = sdpvar(61,10*3);%线路输送功率P

V_angle = sdpvar(54,10*3);

S_pipe_op = sdpvar(22,10*3); %管道运行流量

S_d_g_analog1 = sdpvar(17,10*3);

S_d_g_analog2 = sdpvar(17,10*3);

P_d_analog1 = sdpvar(54,10*3);

P_d_analog2 = sdpvar(54,10*3);

%% ++++++++++++++++++运行子问题++++++++++++++++++

%% ++++++++++++++运行子问题的公共约束条件++++

% +++++++外部接入电源出力约束

F_op_analog = [ P_grid_op <= 15/100*ones(2,10*3)];

F_op_analog = F_op_analog + [ P_grid_op >=0 ];

% +++++++电力潮流约束

for t=1:10

for l=22:61

if B(l,3*t) == 0

F_op_analog = F_op_analog +[ P_line_op(l,3*t-2:3*t) == zeros(1,3) ];

end

% +++++++配电线路流量及容量约束

F_op_analog = F_op_analog + [ V_angle(53,10*3) ==0 ];%平衡节点

F_op_analog = F_op_analog + [ ( P_line_op == B .* ( A_connect_node_line_TR * V_angle ) ):'line' ];

F_op_analog = F_op_analog + [ P_line_op <= S_line_max ];

F_op_analog = F_op_analog + [ P_line_op >= -S_line_max ];

for t=1:10

for i=32:50

if k_node_connect(i,3*t) == 0 %判断节点和线路是否连接，如果不连接，则负荷得不到满足

F_op_analog = F_op_analog + [ V_angle(i,3*t-2:3*t) == zeros(1,3) ] ;

end

% +++++++ 燃气电厂出力约束

F_op_analog = F_op_analog + [ P_DNGFP_op >= 0 ];

F_op_analog = F_op_analog +[ ( P_DNGFP_op <= S_DNGFP_max ) :'DNGFP' ];

% +++++++城市门站出力约束

F_op_analog = F_op_analog + [ S_GasStation_op >= 0];

F_op_analog = F_op_analog + [ ( S_GasStation_op <= S_GasStation_max ): 'GasStation' ];

% +++++++配气主干管道潮流及容量约束

F_op_analog = F_op_analog + [ ( S_pipe_op >= -S_pipe_max ): 'pipelow'];

F_op_analog = F_op_analog + [ ( S_pipe_op <= S_pipe_max ): 'pipeup'];

% +++++++切负荷量约束

for t=1:10

F_op_analog = F_op_analog + [ sum(sum(P_d_cut_op(:,3*t-2:3*t))) <= 0.005 * sum(sum(P_d(:,3*t-2:3*t))) ];

end

F_op_analog = F_op_analog + [ P_d >= P_d_cut_op ];

F_op_analog = F_op_analog + [ 0 <= P_d_cut_op ];

ddd = binvar(54,10*3);

ttt = binvar(54,10*3);

iii = binvar(17,10*3);

jjj = binvar(17,10*3);

for t=1:30

for i=1:54

F_op_analog = F_op_analog + implies(A_connect_node_DNGFP(i,:) * P_DNGFP_op(:,t) + A_connect_node_grid(i,:) * P_grid_op(:,t) + P_d_cut_op(i,t) <= P_d(i,t), ddd(i,t) );

F_op_analog = F_op_analog + implies(ddd(i,t), P_d_analog1(i,t) == A_connect_node_line(i,:) * P_line_op(:,t) + P_d(i,t) - ( A_connect_node_DNGFP(i,:) * P_DNGFP_op(:,t) + A_connect_node_grid(i,:) * P_grid_op(:,t) + P_d_cut_op(i,t) ));

F_op_analog = F_op_analog + implies(ddd(i,t) , P_d_analog2(i,t) == 0 );

F_op_analog = F_op_analog + implies(ddd(i,t) , P_d_analog1(i,t) <= P_d(i,t) - ( A_connect_node_DNGFP(i,:) * P_DNGFP_op(:,t) + A_connect_node_grid(i,:) * P_grid_op(:,t) + P_d_cut_op(i,t) ) );

F_op_analog = F_op_analog + implies(ddd(i,t) , P_d_analog1(i,t) >= 0 );

F_op_analog = F_op_analog + implies(A_connect_node_DNGFP(i,:) * P_DNGFP_op(:,t) + A_connect_node_grid(i,:) * P_grid_op(:,t) + P_d_cut_op(i,t) >= P_d(i,t), ttt(i,t) );

F_op_analog = F_op_analog + implies(ttt(i,t), P_d_analog2(i,t) == -A_connect_node_line(i,:) * P_line_op(:,t) + ( A_connect_node_DNGFP(i,:) * P_DNGFP_op(:,t) + A_connect_node_grid(i,:) * P_grid_op(:,t) + P_d_cut_op(i,t) ) - P_d(i,t));

F_op_analog = F_op_analog + implies(ttt(i,t), P_d_analog1(i,t) == 0 );

F_op_analog = F_op_analog + implies(ttt(i,t), P_d_analog2(i,t) <= ( A_connect_node_DNGFP(i,:) * P_DNGFP_op(:,t) + A_connect_node_grid(i,:) * P_grid_op(:,t) + P_d_cut_op(i,t) ) - P_d(i,t));

F_op_analog = F_op_analog + implies(ttt(i,t), P_d_analog2(i,t) >= 0 );

end

for p=1:17

F_op_analog = F_op_analog + implies(A_connect_gasnode_GasStation(p,:) * S_GasStation_op(:,t) - A_connect_gasnode_DNGFP(p,:) * P_DNGFP_op(:,t)/y_efficiency - S_d_g(p,t)<=0, iii(p,t));

F_op_analog = F_op_analog + implies(iii(p,t), S_d_g_analog1(p,t) + A_connect_gasnode_GasStation(p,:) * S_GasStation_op(:,t)- A_connect_gasnode_DNGFP(p,:) * P_DNGFP_op(:,t)/y_efficiency - S_d_g(p,t) - A_connect_gasnode_pipe(p,:) * S_pipe_op(:,t) == 0 );

F_op_analog = F_op_analog + implies(iii(p,t), S_d_g_analog2(p,t) == 0 );

F_op_analog = F_op_analog + implies(iii(p,t), S_d_g_analog1(p,t) <= - A_connect_gasnode_GasStation(p,:) * S_GasStation_op(:,t) + A_connect_gasnode_DNGFP(p,:) * P_DNGFP_op(:,t)/y_efficiency + S_d_g(p,t));

F_op_analog = F_op_analog + implies(iii(p,t), S_d_g_analog1(p,t) >=0);

F_op_analog = F_op_analog + implies(A_connect_gasnode_GasStation(p,:) * S_GasStation_op(:,t) - A_connect_gasnode_DNGFP(p,:) * P_DNGFP_op(:,t)/y_efficiency - S_d_g(p,t)>= 0, jjj(p,t));

F_op_analog = F_op_analog + implies(jjj(p,t), - S_d_g_analog2(p,t) + A_connect_gasnode_GasStation(p,:) * S_GasStation_op(:,t)- A_connect_gasnode_DNGFP(p,:) * P_DNGFP_op(:,t)/y_efficiency - S_d_g(p,t) - A_connect_gasnode_pipe(p,:) * S_pipe_op(:,t) == 0 );

F_op_analog = F_op_analog + implies(jjj(p,t), S_d_g_analog1(p,t) == 0 );

F_op_analog = F_op_analog + implies(jjj(p,t), S_d_g_analog2(p,t) <= A_connect_gasnode_GasStation(p,:) * S_GasStation_op(:,t) - A_connect_gasnode_DNGFP(p,:) * P_DNGFP_op(:,t)/y_efficiency - S_d_g(p,t));

F_op_analog = F_op_analog + implies(jjj(p,t), S_d_g_analog2(p,t) >=0);

end

Johan Löfberg

unread,

Apr 13, 2016, 6:21:32 AM4/13/16

to YALMIP

Well, I directly note that you use P_d_analog1 is used in an implication, but is not explicitly bounded.

Guori Huang

unread,

Apr 14, 2016, 1:01:36 AM4/14/16

to YALMIP

Thank you! If I use the command of implies, Does it means that I can't get the duals ?

Johan Löfberg

unread,

Apr 14, 2016, 1:17:49 AM4/14/16

to YALMIP

implies introduces binary variables, hence no duals

Guori Huang

unread,

Apr 15, 2016, 2:31:48 AM4/15/16

to YALMIP

could we get the duals from QCP?

Johan Löfberg

unread,

Apr 15, 2016, 3:35:38 AM4/15/16

to YALMIP

Depends on the solver, and how you formulate the QCP. If you write them as simple SOCP cones manually, yes SOCP duals are recovered.

Guori Huang

unread,

Apr 15, 2016, 3:41:34 AM4/15/16

to YALMIP

How about cplex?

Johan Löfberg

unread,

Apr 15, 2016, 4:06:31 AM4/15/16

to YALMIP

Nope. gurobi and mosek works

x = sdpvar(2,1);
C = x'*x <= 100; 
optimize(C,-sum(x),sdpsettings('solver','gurobi','gurobi.qcpdual',1));
dual(C)
optimize(C,-sum(x),sdpsettings('solver','mosek'));
dual(C)
C = cone(x,10);
optimize(C,-sum(x),sdpsettings('solver','gurobi','gurobi.qcpdual',1));
dual(C)
optimize(C,-sum(x),sdpsettings('solver','mosek'));
dual(C)

Guori Huang

unread,

Apr 16, 2016, 1:13:54 AM4/16/16

to YALMIP

Thank you very much. Your advice is very helpful. Another question, could we get the ray by YALMIP?

Johan Löfberg

unread,

Apr 16, 2016, 3:21:20 AM4/16/16

to YALMIP

only if you find it in the numerical output using savesolveroutput

Reply all

Reply to author

Forward

Message has been deleted