Lack of Progress, objective (y-u)'*Q*(y-u)
376 views
Skip to first unread message
Adeleh Mohammadi
Jan 15, 2018, 5:52:42 AM1/15/18
to YALMIP
Dear all,
I have a linear MIMO system with A,B,C,D matrices.
objective of the form: (y-u)'*Q*(y-u)
constraints: range of possible values for one of the inputs and one of the outputs. L < u < U , Ly < y <Uy
Solver:Sdpt3
and I get the verbose that 'Lack of progress'.
Would you please give me your thoughts where the problem could be?
I add that the problem is solved and give results when I remove the constraints.
Can I actually use YALMIP for this form of objective function?
Many thanks
Adeleh
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
num. of constraints = 56
dim. of socp var = 12, num. of socp blk = 1
dim. of linear var = 20
dim. of free var = 45
*** convert ublk to linear blk
********************************************************************************************
SDPT3: homogeneous self-dual path-following algorithms
********************************************************************************************
version predcorr gam expon
HKM 1 0.000 1
it pstep dstep pinfeas dinfeas gap mean(obj) cputime kap tau theta
--------------------------------------------------------------------------------------------
0|0.000|0.000|2.0e+04|3.2e+00|1.8e+05|-2.721954e+05| 0:0:00|1.8e+06|1.0e+00|1.0e+00| chol 1 1
1|0.001|0.001|2.0e+04|3.2e+00|1.8e+05|-2.721503e+05| 0:0:01|1.8e+06|1.0e+00|1.0e+00| chol 1 2
2|0.001|0.001|2.3e+04|3.6e+00|2.1e+05|-3.131912e+05| 0:0:01|1.8e+06|1.0e+00|1.1e+00| chol 1 1
3|0.006|0.006|3.8e+04|5.9e+00|3.7e+05|-5.070290e+05| 0:0:01|1.8e+06|1.0e+00|1.9e+00| chol 1 1
4|0.043|0.043|1.9e+05|2.9e+01|2.5e+06|-2.511153e+06| 0:0:01|1.8e+06|1.0e+00|9.2e+00| chol 1 1
5|1.000|1.000|7.7e+06|1.2e+03|7.5e+08|-1.035134e+08| 0:0:01|1.6e+06|1.0e+00|3.8e+02| chol 1 1
6|0.671|0.671|3.2e+06|5.1e+02|3.1e+08|-4.355323e+07| 0:0:01|5.4e+06|1.0e+00|1.6e+02| chol 1 1
7|0.952|0.952|2.2e+05|3.6e+01|1.9e+07|-2.920561e+06| 0:0:01|6.3e+06|1.0e+00|1.1e+01| chol 1 1
8|1.000|1.000|9.6e+04|1.6e+01|9.2e+06|-1.287988e+06| 0:0:01|1.7e+06|1.0e+00|4.7e+00| chol 1 1
9|1.000|1.000|2.7e+04|4.6e+00|2.5e+06|-3.542189e+05| 0:0:01|3.3e+05|1.0e+00|1.3e+00| chol 1 1
10|0.943|0.943|2.3e+03|3.8e-01|2.0e+05|-2.552636e+04| 0:0:01|7.5e+04|9.0e-01|1.0e-01| chol 1 1
11|0.936|0.936|1.6e+02|2.8e-02|1.3e+04| 5.861599e+02| 0:0:01|9.8e+02|9.7e-01|7.5e-03| chol 1 1
12|0.847|0.847|5.2e+01|1.2e-02|4.2e+03|-1.933814e+02| 0:0:01|1.5e+02|1.1e+00|2.7e-03| chol 1 1
13|1.000|1.000|9.9e+00|8.2e-03|8.5e+02|-6.452446e+01| 0:0:01|3.2e+01|1.1e+00|5.4e-04| chol 1 1
14|0.991|0.991|6.7e-01|7.2e-03|5.1e+01|-6.236824e+00| 0:0:01|7.5e+00|1.1e+00|3.7e-05| chol 1 1
15|0.982|0.982|2.4e-02|6.5e-03|1.4e+00|-3.133262e-01| 0:0:01|6.0e-01|1.1e+00|1.3e-06| chol 1 1
16|0.992|0.992|6.2e-04|5.9e-03|3.2e-02|-9.118175e-03| 0:0:01|2.4e-02|1.1e+00|3.4e-08| chol 1 1
17|1.000|1.000|1.2e-05|5.3e-03|7.0e-04|-1.878572e-04| 0:0:01|5.7e-04|1.1e+00|6.8e-10| chol 1 1
18|0.981|0.981|2.3e-06|2.2e-03|2.1e-04|-5.581393e-05| 0:0:02|3.5e-05|1.1e+00|1.3e-10| chol 1 1
19|0.190|0.190|9.9e-07|1.9e-03|2.1e-04|-1.245541e-04| 0:0:02|2.9e-05|1.1e+00|5.3e-11| chol 1 1
20|0.035|0.035|1.4e-05|1.9e-03|3.1e-04|-5.975876e-04| 0:0:02|2.8e-05|1.1e+00|0.0e+00| chol 1 1
21|0.011|0.011|6.3e-05|1.9e-03|5.9e-04|-2.125171e-03| 0:0:02|2.8e-05|1.1e+00|0.0e+00|
Stop: progress is too slow
-------------------------------------------------------------------
number of iterations = 21
primal objective value = -4.24839012e-03
dual objective value = -1.95274294e-06
gap := trace(XZ) = 5.93e-04
relative gap = 5.91e-04
actual relative gap = -4.23e-03
rel. primal infeas = 6.31e-05
rel. dual infeas = 1.87e-03
norm(X), norm(y), norm(Z) = 7.1e-01, 3.7e+07, 3.2e+04
norm(A), norm(b), norm(C) = 8.6e+07, 1.0e+00, 6.1e+07
Total CPU time (secs) = 1.59
CPU time per iteration = 0.08
termination code = -5
DIMACS: 6.3e-05 0.0e+00 1.9e-03 0.0e+00 -4.2e-03 5.9e-04
-------------------------------------------------------------------
ans =
'Lack of progress '
Erling D. Andersen
Jan 15, 2018, 6:46:27 AM1/15/18
to YALMIP
I can see norm of C and A are very large so your model is likely to be badly scaled.
You may also want to try MOSEK. If that fails too, then it is very likely your problem that is broken.
Johan Löfberg
Jan 15, 2018, 6:47:12 AM1/15/18
to YALMIP
Why are you using SDPT3 to solve a simple quadratic program?
You should install a dedicated QP solver (or use it if you already have one installed)
Erling D. Andersen
Jan 15, 2018, 6:56:50 AM1/15/18
to YALMIP
Using SOCP solver is not that stupid. From a theoretical point of view you cannot do any better computational complexity wise. Please correct me if I wrong about this.
Also it is very common among the commercial QP solvers to convert it to conic form internally and solve it that way.
Johan Löfberg
Jan 15, 2018, 7:12:31 AM1/15/18
to YALMIP
True, its about implementation, what I meant is more that I've seen to many poorly solved QPs in sdpt3/sedumi (+ you have the active set option when you stay in the QP world)
Adeleh Mohammadi
Jan 16, 2018, 4:43:33 AM1/16/18
to YALMIP
Thank you for your reply Erling.
I used another solver and it gives me the same issue, so as you suggested it's my problem formulation that is broken. It is a control of a chiller, I put the code here, what do you think about the problem formulation?
I only have some of the inputs and outputs involved in the cost and constraints, so I considered the rest as disturbance and didn't include them among the decision variables.
%% MPC of chiller
clc
close all
clear all
yalmip('clear')
%% Model Parameters and State-space matrices
data = xlsread('09012018data.csv');
% sampling data every 5 minutes
NoM = 5;
sampling_period = NoM*60;
sampled_data = [];
for i = 1:sampling_period:length(data)
sampled_data = [sampled_data;data(i,:)];
end
Time = sampled_data(:,1);
dTsetE = sampled_data(:,2);
dFlow = sampled_data(:,3);
dTiE = sampled_data(:,4);
dTiC = sampled_data(:,5);
dToutE = sampled_data(:,6);
dPout = sampled_data(:,7);
dCOP = sampled_data(:,8);
A = [-0.055643256 -2.79E-07 0.000971282 1 0 0
10354.94655 0.223207545 356.7890554 0 1 0
-6.467478387 2.81E-05 0.022706046 0 0 1
-0.014415444 2.60E-07 0.002964755 0 0 0
-812.961162 0.022406063 -30.49280997 0 0 0
0.018877885 -5.01E-06 0.084396962 0 0 0];
B = [0.017021927 -0.001889088 -0.015300244 -0.007369483
2226.598087 -190.7771176 -348.0510377 -89.55780513
-0.030880581 -0.009586317 0.286861591 -0.025348384
-0.003255077 -0.001385605 -0.012293342 7.07E-05
827.9699533 -168.1696995 -1583.890508 -6.360500619
-0.531436147 0.050092133 0.495164075 0.004485626];
B1 = [0.017021927 -0.001889088
2226.598087 -190.7771176
-0.030880581 -0.009586317
-0.003255077 -0.001385605
827.9699533 -168.1696995
-0.531436147 0.050092133]; % corresponding to Tset and Flow
B2 = [-0.015300244 -0.007369483
-348.0510377 -89.55780513
0.286861591 -0.025348384
-0.012293342 7.07E-05
-1583.890508 -6.360500619
0.495164075 0.004485626]; % corresponding to TiE, ToC
C = [1 0 0 0 0 0;
0 1 0 0 0 0;
0 0 1 0 0 0];
Cr = [1 0 0 0 0 0;
0 1 0 0 0 0];
D = [0.973267178 0.016742346 0.035864605 0.006784957
68966.72868 -6630.599254 -66807.43608 35.43101393
-2.424358679 0.19654465 1.896918069 0.053732785];
Dr = [0.973267178 0.016742346 0.035864605 0.006784957
68966.72868 -6630.599254 -66807.43608 35.43101393];
Dr1 = [0.973267178 0.016742346
68966.72868 -6630.599254];
Dr2 = [ 0.035864605 0.006784957
-66807.43608 35.43101393];
D1 = [0.973267178 0.016742346
68966.72868 -6630.599254
-2.424358679 0.19654465 ];
D2 = [ 0.035864605 0.006784957
-66807.43608 35.43101393
1.896918069 0.053732785];
nx = length(A); % Number of states
[ny , nu] = size(D); % Number of inputs and outputs
%% MPC problem
N = 5; % MPC horizon
% output objective 'ToutEvaporator','Power','COP'
Wq = 1;
Q = Wq*eye(ny); % Input weighting matrix
% Output reference
ToEref = 281.15; % Set to be 8 C
Pref = 0; % P to reach as close as possible to 0
% Output Constraints
ToEmax = 285.15; % This is the max value in physical system (12 C)
ToEmin = 275.15; % This is the min value in physical system (2 C)
% Pmax = max(dPout); % max value of data
Pmax = 10e6;
% Pmin = min(dPout); % min value of data
Pmin = 0; % min value of data
% Input objective weighting
Wr = 1; % state weighting gain
R = Wr*eye(nu); % state weighting matrix
% Input constraints
TsetEmax = max(dTsetE);
TsetEmin = min(dTsetE);
Flowmax = max(dFlow);
% Flowmin = min(dFlow);
Flowmin = 0;
%% Optimization problem: Decision variables, objective, constraints
% Vector of inputs including decision variables
ut = sdpvar(repmat(2,1,N),repmat(1,1,N)); % TsetE and Flow
x = sdpvar(repmat(nx,1,N+1),repmat(1,1,N+1)); % states
yt = sdpvar(repmat(2,1,N+1),repmat(1,1,N+1)); % outputs
constraints = [];
objective = 0;
for h = 1:N
% disturbance vector [Load;weather]
w = [dTiE(h);dTiC(h)];
% Minimization of the term (TsetE - ToE)
objective = objective + (ut{h}(1,:)-yt{h}(1,:))'*1*(ut{h}(1,:)-yt{h}(1,:));
% Power minimization
objective = objective + (yt{h}(2,:)-Pref)'*1*(yt{h}(2,:)-Pref);
% Minimization of Tset, mFlow
% objective = objective + (ut{h}(1,:) - ur)'*R*(ut{h}(1,:) - ur);
% System Dynamics
constraints = [constraints, x{h+1} == A*x{h} + B1*ut{h} + B2*w];
constraints = [constraints, yt{h} == Cr*x{h} + Dr1*ut{h} + Dr2*w];
constraints = [constraints, TsetEmin<= ut{h}(1,:) <= TsetEmax];
constraints = [constraints, Flowmin<= ut{h}(2,:) <= Flowmax];
constraints = [constraints, Pmin <= yt{h}(2,:) <= Pmax];
end
ops = sdpsettings('verbose',2,'solver','sdpt3'); % 'qpoases' , 'mosek'
controller = optimizer(constraints,objective,ops,x{1},[ut{:}]);
x = zeros(nx,1);
control_history = [];
state_history = [];
output_history = [];
Nsim = length(sampled_data);
for i = 1:Nsim
[U,errorcode] = controller{x};
if errorcode
yalmiperror(errorcode)
break
end
w = [dTiE(h);dTiC(h)];
% y = [yt{h};dCOP(h)];
Uf = [U(:,1);w];
x = A*x + B1*U(:,1) + B2*w;
y = C*x + D1*U(:,1) + D2*w;
control_history = [control_history,Uf];
state_history = [state_history,x];
output_history = [output_history,y];
end
value(objective)
Adeleh Mohammadi
Jan 16, 2018, 4:46:55 AM1/16/18
to YALMIP
Dear Johan,
Thank you very much for your comments and thoughts.
Would you please give me your thoughts on the problem formulation in the full code that I just posted? This is an MPC for a chiller.
Thank you in advance,
Adeleh
Adeleh Mohammadi
Jan 16, 2018, 4:49:13 AM1/16/18
to YALMIP
and the results when I run:
num. of constraints = 51
dim. of socp var = 12, num. of socp blk = 1
dim. of linear var = 30
dim. of free var = 40
*** convert ublk to linear blk
********************************************************************************************
SDPT3: homogeneous self-dual path-following algorithms
********************************************************************************************
version predcorr gam expon
HKM 1 0.000 1
it pstep dstep pinfeas dinfeas gap mean(obj) cputime kap tau theta
--------------------------------------------------------------------------------------------
0|0.000|0.000|2.3e+04|6.7e+00|7.3e+09| 2.819481e+08| 0:0:00|7.3e+10|1.0e+00|1.0e+00| chol 1 2
1|0.000|0.000|2.3e+04|6.7e+00|7.3e+09| 2.817833e+08| 0:0:01|7.3e+10|1.0e+00|1.0e+00| chol 2 3
2|0.000|0.000|4.4e+04|1.3e+01|7.8e+09| 1.447154e+08| 0:0:01|6.8e+10|1.0e+00|1.9e+00|
Stop: steps are too short
-------------------------------------------------------------------
number of iterations = 2
primal objective value = 6.11782990e+08
dual objective value = 0.00000000e+00
gap := trace(XZ) = 8.65e+09
relative gap = 2.83e+01
actual relative gap = 1.00e+00
rel. primal infeas = 2.50e+04
rel. dual infeas = 7.21e+00
norm(X), norm(y), norm(Z) = 6.5e+01, 0.0e+00, 1.3e+08
norm(A), norm(b), norm(C) = 1.5e+08, 1.0e+00, 6.5e+07
Total CPU time (secs) = 1.28
CPU time per iteration = 0.64
termination code = -5
DIMACS: 2.5e+04 0.0e+00 7.2e+00 0.0e+00 1.0e+00 1.4e+01
-------------------------------------------------------------------
ans =
'Lack of progress '
ans =
NaN
Erling D. Andersen
Jan 16, 2018, 4:57:48 AM1/16/18
to YALMIP
I have comments on your model but I suggest you try out MOSEK from https://mosek.com.
Adeleh Mohammadi
Jan 16, 2018, 5:05:14 AM1/16/18
to YALMIP
Dear Erling,
I downloaded MOSEK a while ago and got the academic license with my student email, I installed it on my Windows device and placed the license file in the home directory as suggested,
This means on Windows the license file should normally be save to:
c:\users\<userid>\mosek\mosek.lic
c:\users\<userid>\mosek\mosek.lic
However when I use MOSEK as the solver
ops = sdpsettings('verbose',2,'solver','mosek');
it says:
Warning: Solver not found (mosek)
Error using optimizer (line 193)
Failed exporting the model: Solver not found (mosek)
Error in toyalmipdevelopers (line 152)
controller = optimizer(constraints,objective,ops,x{1},[ut{:}]);
How could I resolve this?
Erling D. Andersen
Jan 16, 2018, 5:18:05 AM1/16/18
to YALMIP
The last stuff has to do with yalmip and that most likely means you did not follow the instructions. So let us start by get mosek running as a stand alone.
and followed the instructions.
What happens when you run
mosekdiag
?
Adeleh Mohammadi
Jan 16, 2018, 5:28:53 AM1/16/18
to YALMIP
aha, now it works, I used to add the main path of MOSEK rather than the MATLAB folder.
it says:
Matlab version: 9.3.0.713579 (R2017b)
Architecture : PCWIN64
The mosek optimizer executed successfully from the command line:
MOSEK Version 8.1.0.34 (Build date: 2017-11-28 12:38:49)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86
FlexLM
Version : 11.14
Hostname : IRUSE_LAPTOP
Host ID : "00ff739ce812 9cb70d5505dd 0aedb9295df7 d067e54ba37e 08edb9295df7"
License path : C:\Users\IRUSE\mosek\mosek.lic
Operating system variables
PATH :
*** No input file specfied. No optimization is performed.
Return code - 0 [MSK_RES_OK]
mosekopt: D:\Users\IRUSE\Desktop\Solvers\mosek\8\toolbox\r2014aom\mosekopt.mexw64
mosekopt is working correctly.
Warning: MOSEK Fusion is not configured correctly; check that mosek.jar is added to the javaclasspath.
> In mosekdiag (line 72)
and when I run the same code, it says:
Optimizer - threads : 2
Optimizer - solved problem : the primal
Optimizer - Constraints : 26
Optimizer - Cones : 2
Optimizer - Scalar variables : 58 conic : 28
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 344 after factor : 351
Factor - dense dim. : 0 flops : 8.22e+03
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.0e+00 3.1e+05 1.0e+00 0.00e+00 1.280000000e+02 0.000000000e+00 1.0e+00 0.13
1 2.5e-01 7.9e+04 1.3e-01 -1.00e+00 -4.832919573e+04 -3.525410343e+02 2.5e-01 0.33
2 6.1e-02 1.9e+04 1.5e-02 -1.00e+00 -2.533467325e+05 -3.027306409e+03 6.1e-02 0.34
3 1.8e-02 5.7e+03 2.5e-03 -1.00e+00 -8.907405246e+05 -1.367175672e+04 1.8e-02 0.36
4 6.7e-03 2.1e+03 5.5e-04 -9.98e-01 -2.488524169e+06 -5.567730304e+04 6.7e-03 0.36
5 2.0e-03 6.3e+02 9.2e-05 -9.96e-01 -8.237968943e+06 -2.519869259e+05 2.0e-03 0.38
6 6.0e-04 1.9e+02 1.5e-05 -9.87e-01 -2.719504694e+07 -1.016968925e+06 6.0e-04 0.39
7 2.1e-04 6.4e+01 3.1e-06 -9.61e-01 -7.550009562e+07 -3.395641406e+06 2.1e-04 0.39
8 8.9e-05 2.8e+01 9.4e-07 -9.03e-01 -1.585295721e+08 -1.041523813e+07 8.9e-05 0.41
9 5.2e-05 1.6e+01 4.4e-07 -8.51e-01 -2.500334482e+08 -1.285074687e+07 5.2e-05 0.42
10 1.6e-05 5.0e+00 9.1e-08 -7.97e-01 -5.994195119e+08 -8.914510112e+07 1.6e-05 0.42
11 9.7e-06 3.0e+00 4.8e-08 -4.99e-01 -8.149122070e+08 -1.185063993e+08 9.7e-06 0.44
12 1.5e-06 4.8e-01 4.4e-09 -5.65e-01 -2.394690887e+09 -4.161034959e+08 1.5e-06 0.44
13 4.4e-07 1.4e-01 9.6e-10 -4.31e-01 -4.314810628e+09 -9.223882680e+08 4.4e-07 0.45
14 1.3e-07 3.9e-02 3.5e-10 -8.48e-02 -4.500471779e+09 -2.308892177e+09 1.3e-07 0.47
15 3.6e-08 1.1e-02 6.4e-11 -2.67e-01 -9.466577249e+09 -4.186993985e+09 3.6e-08 0.48
16 5.3e-09 1.7e-03 2.7e-11 2.34e-01 -1.089514557e+10 -1.022745889e+10 5.3e-09 0.48
17 1.8e-09 5.6e-04 4.9e-12 -2.92e-01 -1.681358621e+10 -1.458652838e+10 1.8e-09 0.50
18 8.1e-10 2.1e-04 6.1e-12 9.95e-01 -2.132561846e+10 -2.112560977e+10 6.7e-10 0.52
19 9.0e-11 2.3e-05 3.4e-12 1.33e+00 -1.984739319e+10 -1.983867666e+10 7.3e-11 0.53
20 4.4e-11 1.6e-06 1.2e-12 1.10e+00 -1.929503219e+10 -1.929448363e+10 4.9e-12 0.53
21 4.4e-11 1.6e-06 3.4e-12 1.00e+00 -1.929502121e+10 -1.929447307e+10 4.9e-12 0.55
22 2.5e-11 2.5e-07 2.4e-11 1.00e+00 -1.928360439e+10 -1.928349211e+10 8.0e-13 0.56
23 7.4e-12 4.1e-08 3.5e-11 1.00e+00 -1.928511052e+10 -1.928507815e+10 1.3e-13 0.56
24 3.6e-12 5.5e-08 3.6e-10 1.00e+00 -1.928495669e+10 -1.928494373e+10 4.0e-14 0.58
25 2.3e-12 1.3e-07 8.3e-10 1.00e+00 -1.928472489e+10 -1.928471807e+10 1.5e-14 0.59
26 9.8e-13 2.0e-07 2.2e-09 1.00e+00 -1.928477377e+10 -1.928477081e+10 5.5e-15 0.59
27 9.8e-13 2.0e-07 2.2e-09 1.00e+00 -1.928477377e+10 -1.928477081e+10 5.5e-15 0.61
Optimizer terminated. Time: 0.73
Interior-point solution summary
Problem status : UNKNOWN
Solution status : UNKNOWN
Primal. obj: -1.9284773765e+10 nrm: 2e+10 Viol. con: 1e+01 var: 0e+00 cones: 0e+00
Dual. obj: -1.9284770810e+10 nrm: 2e+10 Viol. con: 0e+00 var: 2e+02 cones: 0e+00
Optimizer summary
Optimizer - time: 0.73
Interior-point - iterations : 28 time: 0.63
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00
ans =
'Unknown problem in solver (Turn on 'debug' in sdpsettings) '
Erling D. Andersen
Jan 16, 2018, 6:21:23 AM1/16/18
to YALMIP
The optimal primal and dual solutions are about of the size 10^10 i.e. ||x|| is 10^10.
So you should ask yourself why is the optimal primal and dual solutions very large.
Reply all
Reply to author
Forward
0 new messages