Suboptimal termination and numerical issues
Answered
For the following code in Matlab, the solver runs into numerical trouble and terminates sub-optimally. I referred the Matrix coefficient range described on https://www.gurobi.com/documentation/9.5/refman/does_my_model_have_numeric.html but this does not appear to be the problem. I will appreciate any advisable steps. The code is below:
% boundary conditions
boundary.v0 = 1; % [m/s] initial velocity
boundary.s0 = 0; % [m] initial position
boundary.E0 = 306e6; % [J] initial battery stored energy
boundary.v1 = 1; % [m/s] velocity at final time
boundary.s1 = 6000; % [m] distance to travel
boundary.tf = 200; % [s] travel time
% vehicle parameters
vehiclepar.cr = 0.01; % [N/kN] rolling resistance
vehiclepar.cd = 0.24; % [-] aerodynamic drag coefficient
vehiclepar.A = 2.34; % [m2] frontal area
vehiclepar.rho = 1.2; % [kg/m3] air density
vehiclepar.g = 9.81; % [m/s2] gravitational constant
vehiclepar.m = 2108; % [kg] vehicle mass
vehiclepar.Pd = -20000; % [W] brake power
vehiclepar.Pm_up = 310e3; % [W] max driving power of EM
vehiclepar.Pm_lw = -vehiclepar.Pm_up; % [W] max braking power of EM
vehiclepar.gm_m = 0.9; % [-] electric machine efficiency
vehiclepar.R = 0.32; % [Ohm] battery internal resistance
vehiclepar.U0 = 346; % [V] discharged battery voltage
vehiclepar.Ebat = 306e6; % [J] battery capacity
vehiclepar.roadgrad= 0; % [deg] road gradient
yalmip('clear')
delta_s = 2; % s1/Δs should be an integer
Ek = sdpvar(boundary.s1/delta_s+1,1); % state variable 1 of 3
Es = sdpvar(boundary.s1/delta_s+1,1); % state variable 2 of 3
t = sdpvar(boundary.s1/delta_s+1,1); % state variable 3 of 3
v = sdpvar(boundary.s1/delta_s+1,1);
delta_t = sdpvar(boundary.s1/delta_s,1); % Δt
Fs = sdpvar(boundary.s1/delta_s,1); % force equivalent of battery power before internal losses
Fb = sdpvar(boundary.s1/delta_s,1); % force equivalent of battery power after internal losses
Fl = sdpvar(boundary.s1/delta_s,1); % rolling resistance + aero drag
Fm = sdpvar(boundary.s1/delta_s,1); % force equivalent of propulsive power motor provides
Fd = sdpvar(boundary.s1/delta_s,1); % brake force
objfun = Fs'*delta_s*ones(boundary.s1/delta_s,1);
cons = [ v(1) == boundary.v0; % initial vel
v(end) == boundary.v1; % final vel
t(1) == 0; % initial time
sum(delta_t) == boundary.tf; % final time
Ek >= 1/2*vehiclepar.m*v.^2; % Ek=0.5*m*v^2
cone([ delta_t/delta_s+v(2:end);
2*ones(size(delta_t));
delta_t/delta_s-v(2:end) ]); % v*Δt/Δs = 1
Fb >= Fm./vehiclepar.gm_m; % Fb>Fm/η
Fb >= Fm.*vehiclepar.gm_m; % Fb>Fm*η
cone([ Fs-Fb+delta_t/delta_s;
2*sqrt(vehiclepar.R)/vehiclepar.U0*Fs;
Fs-Fb+delta_t/delta_s ]) ]; % (Fs-Fb)*Δt/Δs>=Fs^2*R/U0^2
for i = 1:boundary.s1/delta_s
cons = [cons;
Fl(i) == vehiclepar.cr*vehiclepar.m*vehiclepar.g ...
+ vehiclepar.rho*vehiclepar.A*vehiclepar.cd/vehiclepar.m*Ek(i); % rolling resistance + aero drag
Ek(i+1) == Ek(i) + (Fm(i) - Fd(i) - Fl(i))*delta_s; % Ek state constraint
Es(i+1) == Es(i) - Fs(i)*delta_s; % Es state constraint
t(i+1) == t(i) + delta_t(i) ]; % time state constraint
end
opts = sdpsettings( 'debug',1, ...
'solver','gurobi', ...
'gurobi.NonConvex',0, ... % error if quadratic constraints are non-convex
'gurobi.Method',2, ... % barrier for QCP
'gurobi.NumericFocus',3, ... % higher preference for accuracy over speed
'showprogress',1 );
sol = optimize(cons,objfun,opts)
-
The residuals increase in the output:
+ Solver chosen : GUROBI-GUROBI+ Processing objective function+ Processing constraints+ Calling GUROBI
Academic license - for non-commercial use only - expires 2022-06-18
Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (win64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 45007 rows, 57007 columns and 123009 nonzeros
Model fingerprint: 0xdec0ed1b
Model has 3003 quadratic constraints
Coefficient statistics:
Matrix range [3e-04, 3e+01]
QMatrix range [1e+00, 1e+00]
Objective range [2e+00, 2e+00]
Bounds range [0e+00, 0e+00]
RHS range [5e-01, 2e+02]
Presolve removed 21005 rows and 18006 columns
Presolve time: 0.13s
Presolved: 24005 rows, 39004 columns, 69003 nonzeros
Presolved model has 3003 second-order cone constraints
Ordering time: 0.00s
Barrier statistics:
Dense cols : 2
Free vars : 11999
AA' NZ : 6.600e+04
Factor NZ : 1.800e+05 (roughly 27 MBytes of memory)
Factor Ops : 1.477e+06 (less than 1 second per iteration)
Threads : 1
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 -3.99998400e-04 0.00000000e+00 1.01e+03 1.00e+03 2.22e+05 0s
1 -4.54525687e+08 -6.42426914e+05 5.99e+02 5.95e+02 1.37e+05 0s
2 -8.64227254e+08 -1.16649288e+06 2.27e+02 2.09e+02 6.42e+04 1s
3 -9.32084500e+08 -1.43906621e+06 5.82e+01 4.95e+01 4.65e+04 1s
4 -1.16068600e+09 -1.44821872e+06 2.35e+01 1.69e+01 2.17e+04 1s
5 -1.28103857e+09 -1.41411423e+06 1.19e+01 9.75e+00 9.94e+03 1s
6 -8.58073869e+08 -1.31300471e+06 4.82e+00 4.79e+00 4.53e+03 1s
7 -6.18524094e+08 -1.19188077e+06 2.36e+00 6.59e+00 3.44e+03 1s
8 -5.57858882e+08 -1.08281864e+06 3.38e+00 1.98e+01 2.29e+03 1s
9 -9.70215990e+07 -8.14599794e+05 9.32e-01 6.97e+00 4.08e+03 1s
10 -9.22608332e+07 -6.91108363e+05 7.82e-01 1.77e+01 2.94e+03 2s
11 -3.99798410e-04 0.00000000e+00 5.14e+03 5.00e+03 5.56e+06 2s
12 -1.02327024e+08 -5.09747690e+07 4.64e+03 4.65e+03 5.36e+07 2s
13 -4.86431011e+08 -3.36849295e+08 3.43e+03 3.38e+03 2.19e+08 2s
14 -9.22062502e+08 -7.57862052e+08 1.55e+03 1.45e+03 2.07e+08 2s
15* -1.03920395e+09 -9.50532013e+08 3.30e+02 4.41e+02 8.70e+03 2s
16* -1.11880042e+09 -1.08003111e+09 1.09e+02 1.47e+02 2.92e+03 2s
17* -1.16221797e+09 -1.15455630e+09 2.77e+01 3.41e+01 7.07e+02 3s
18* -1.25638855e+09 -1.25202495e+09 7.81e+00 9.41e+00 1.98e+02 3s
19* -1.26675505e+09 -1.26408016e+09 3.43e+00 4.27e+00 8.94e+01 3s
20 -2.96321607e+08 -4.62898238e+08 1.33e+00 4.38e+00 1.79e+05 3s
21 -2.46864449e+08 -3.11774023e+08 2.35e+00 8.17e+00 9.31e+04 3s
22* -2.19273371e+08 -2.46462529e+08 1.99e-01 4.39e-01 4.90e+00 3s
23* -5.53638836e+07 -1.36452008e+08 8.43e-02 3.61e+00 1.22e+00 4s
24* -9.42201660e+07 -1.36390751e+08 3.89e-02 2.24e+00 2.60e-01 4s
25* -1.09042715e+08 -1.40895312e+08 7.55e-03 5.93e+00 1.92e-01 4s
26* -7.59251996e+07 -1.05668911e+08 2.00e-01 1.66e+01 5.47e-02 4s
27* -8.05097206e+07 -9.14561050e+07 1.29e-01 7.74e+01 1.23e-02 4s
28* -1.17526437e+08 -1.19701577e+08 1.18e+00 1.75e+02 1.79e-03 4s
29* -1.14285395e+08 -1.16336882e+08 4.70e+00 9.89e+01 8.35e-04 5s
30* -1.27795421e+08 -1.16355808e+08 3.87e+02 3.17e+01 3.14e-04 5s
31* -1.27795435e+08 -1.16355843e+08 3.87e+02 3.17e+01 3.13e-04 5s
32* -1.27795458e+08 -1.16355866e+08 3.87e+02 3.17e+01 3.12e-04 5s
33* -1.27795465e+08 -1.16355895e+08 3.87e+02 3.17e+01 3.12e-04 5s
34* -1.27795471e+08 -1.16355942e+08 3.87e+02 3.17e+01 3.12e-04 5s
35* -1.27795493e+08 -1.16355986e+08 3.87e+02 3.17e+01 3.11e-04 5s
36* -1.27795498e+08 -1.16356014e+08 3.87e+02 3.17e+01 3.11e-04 5s
37* -1.27795499e+08 -1.16356017e+08 3.87e+02 3.17e+01 3.11e-04 6s
38* -1.27795500e+08 -1.16356030e+08 3.87e+02 3.17e+01 3.11e-04 6s
39* -1.27795504e+08 -1.16356043e+08 3.87e+02 3.17e+01 3.11e-04 6s
40* -1.27795504e+08 -1.16356053e+08 3.87e+02 3.17e+01 3.11e-04 6s
41* -1.27795507e+08 -1.16356061e+08 3.87e+02 3.17e+01 3.10e-04 6s
42* -1.27795508e+08 -1.16356070e+08 3.87e+02 3.17e+01 3.10e-04 6s
43* -1.27795508e+08 -1.16356082e+08 3.87e+02 3.17e+01 3.11e-04 6s
44* -1.27795509e+08 -1.16356086e+08 3.87e+02 3.17e+01 3.10e-04 6s
45* -1.27795509e+08 -1.16356096e+08 3.87e+02 3.17e+01 3.10e-04 7s
46* -1.27795511e+08 -1.16356102e+08 3.87e+02 3.17e+01 3.10e-04 7s
47* -1.27795514e+08 -1.16356114e+08 3.87e+02 3.17e+01 3.10e-04 7s
48* -1.27795514e+08 -1.16356123e+08 3.87e+02 3.17e+01 3.09e-04 7s
49* -1.27795514e+08 -1.16356128e+08 3.87e+02 3.17e+01 3.09e-04 7s
50* -1.27795520e+08 -1.16356145e+08 3.87e+02 3.17e+01 3.09e-04 7s
Barrier performed 50 iterations in 7.26 seconds
Sub-optimal termination - objective -1.27795520e+080 -
Could you try updating to Gurobi's latest version 9.5.x?
Could you try using the BarHomogeneous algorithm?
Could you share an LP or MPS file of your model? The Knowledge Base article How do I write an MPS model file from a third-party API? describes a way of obtaining a model file. To share a model file in the community please follow the instructions in Posting to the Community Forum.
0 -
I updated to the latest version and implemented BarHomogeneous algorithm. Unfortunately, the problem still persists. Please find the MPS model here: https://www.filemail.com/d/apmkkjyeakmiirv
For reference, I have added the code too:
%%System Parameters
% boundary conditions
boundary.v0 = 1; % [m/s] initial velocity
boundary.s0 = 0; % [m] initial position
boundary.E0 = 306e6; % [J] initial battery stored energy
boundary.v1 = 1; % [m/s] velocity at final time
boundary.s1 = 6000; % [m] distance to travel
boundary.tf = 200; % [s] travel
% vehicle parameters
vehiclepar.cr = 0.01; % [N/kN] rolling resistance
vehiclepar.cd = 0.24; % [-] aerodynamic drag coefficient
vehiclepar.A = 2.34; % [m2] frontal area
vehiclepar.rho = 1.2; % [kg/m3] air density
vehiclepar.g = 9.81; % [m/s2] gravitational constant
vehiclepar.m = 2108; % [kg] vehicle mass
vehiclepar.Pd = -20000; % [W] brake power
vehiclepar.Pm_up = 310e3; % [W] max driving power of EM
vehiclepar.Pm_lw = -vehiclepar.Pm_up; % [W] max braking power of EM
vehiclepar.gm_m = 0.9; % [-] electric machine efficiency
vehiclepar.R = 0.32; % [Ohm] battery internal resistance
vehiclepar.U0 = 346; % [V] discharged battery voltage
vehiclepar.Ebat = 306e6; % [J] battery capacity
vehiclepar.roadgrad = 0; % [deg] road gradient
yalmip('clear')
delta_s = 2; % s1/Δs should be an integer
Ek = sdpvar(boundary.s1/delta_s+1,1); % state variable 1 of 3
Es = sdpvar(boundary.s1/delta_s+1,1); % state variable 2 of 3
t = sdpvar(boundary.s1/delta_s+1,1); % state variable 3 of 3
v = sdpvar(boundary.s1/delta_s+1,1);
delta_t = sdpvar(boundary.s1/delta_s,1); % Δt
Fs = sdpvar(boundary.s1/delta_s,1); % force equivalent of battery power before internal losses
Fb = sdpvar(boundary.s1/delta_s,1); % force equivalent of battery power after internal losses
Fl = sdpvar(boundary.s1/delta_s,1); % rolling resistance + aero drag
Fm = sdpvar(boundary.s1/delta_s,1); % force equivalent of propulsive power motor provides
Fd = sdpvar(boundary.s1/delta_s,1); % brake force
objfun = Fs'*delta_s*ones(boundary.s1/delta_s,1);
cons = [ Ek >= 0; % bounds
t >= 0; % bounds
v >= 0; % bounds
Fl>= 0; % bounds
Fd >= 0; % bounds
v(1) == boundary.v0; % initial vel
v(end) == boundary.v1; % final vel
t(1) == 0; % initial time
sum(delta_t) == boundary.tf; % final time
Ek >= 1/2*vehiclepar.m*v.^2; % Ek=0.5*m*v^2
cone([ delta_t/delta_s+v(2:end);
2*ones(size(delta_t));
delta_t/delta_s-v(2:end) ]); % v*Δt/Δs = 1
Fb >= Fm./vehiclepar.gm_m; % Fb>Fm/η
Fb >= Fm.*vehiclepar.gm_m; % Fb>Fm*η
cone([ Fs-Fb+delta_t/delta_s;
2*sqrt(vehiclepar.R)/vehiclepar.U0*Fs;
Fs-Fb+delta_t/delta_s ]) ]; % (Fs-Fb)*Δt/Δs>=Fs^2*R/U0^2
for i = 1:boundary.s1/delta_s
cons = [cons;
Fl(i) == vehiclepar.cr*vehiclepar.m*vehiclepar.g ...
+ vehiclepar.rho*vehiclepar.A*vehiclepar.cd/vehiclepar.m*Ek(i); % rolling resistance + aero drag
Ek(i+1) == Ek(i) + (Fm(i) - Fd(i) - Fl(i))*delta_s; % Ek state constraint
Es(i+1) == Es(i) - Fs(i)*delta_s; % Es state constraint
t(i+1) == t(i) + delta_t(i) ]; % time state constraint
end
opts = sdpsettings( 'debug',1, ...
'solver','gurobi', ...
'gurobi.Method',2, ... % barrier for QCP
'gurobi.NumericFocus',3, ... % higher preference for accuracy over speed
'gurobi.BarHomogeneous',1, ... % useful for recognizing infeasibility or unboundedness
'showprogress',1 );
sol = optimize(cons,objfun,opts)0 -
Update. I tried with the following settings too but they did not solve the problem unfortunately. Do you advise me to try anything else?
opts = sdpsettings( 'debug',1, ...
'solver','gurobi', ... % suitable for second order conic programming
'gurobi.Presolve',0, ... % 0=off, 1=conservative, 2=aggressive
'gurobi.Method',2, ... % barrier for QCP
'gurobi.NumericFocus',3, ... % higher preference for accuracy over speed
'gurobi.ScaleFlag',3, ... % model is scaled to improve the numerical properties of the constraint matrix
'gurobi.CrossoverBasis',1, ... % produces a more numerically stable start basis
'gurobi.BarHomogeneous',1, ... % useful for recognizing infeasibility or unboundedness
'gurobi.FeasibilityTol',0.01, ... % all constraints are satisfied to a tolerance of FeasibilityTol
'gurobi.OptimalityTol',0.01, ... % model will be declared optimal under this value
'showprogress',1 );0 -
I don't think that setting some parameters is enough in this case. I have the suspicion that it is a numerical issue but I am still investigating which may take some time.
In the meantime, could you provide a mathematical formulation of the model? This means in terms of written equations, e.g., in PDF format, instead of computer code. This might make analyzing the issue easier.
0 -
Please find the mathematical formulation of the model at https://www.filemail.com/d/ibaxjcixhyqamyx
I tried to normalise the variables but that deteriorated the coefficient ranges. Let me know the changes you suggest me to make to the model.
0 -
Hi Tushar,
Could you please try the following reformulation and share the mps/lp file here?
Use equality constraints
R0: - 0.5 C0 + C30004 = 0.5
R1: - 32.46536616149585 C9003 + C30005 = 0
R2: 0.5 C0 + C30006 = 0.5to substitute C30004-C30006 and reformulate
qc0: [ - C30004 ^2 + C30005 ^2 + C30006 ^2 ] <= 0
into
qc0: - C0 + [ 1054 C9003 ^2 ] <= 0
This can be done for all quadratic constraints qc0 - qc3000.
Then, use equality constraints R9003-R18002 to reformulate qc3001 by substituting all auxiliary variables C39007-C48006.
Finally, use quality constraints R18003-R27002 to reformulate qc3002 by substituting all auxiliary variables C48007-C57006.
This might better give a better indication where the numerical issues are coming from.
Best regards,
Jaromił0 -
Could you elaborate on what reformulation we are doing here? I do not understand the terms. Are C0, C30004, C9003, C30005, C30006 and C30004 are the new variables that I need to define and add the above constraints into my problem?
0 -
Variables C0, C30004, ... are variables from the MPS file you provided in LP format. You can generate an LP file out of your mps file via
./gurobi_cl resultfile=myLP.lp gurobi.mps
You will then find the above variables and constraints in the myLP.lp file.
In the reformulation I am substituting
C30004 = 0.5 +0.5 C0
C30005 = 32.46536616149585 C9003
C30006 = 0.5 - 0.5 C0into
qc0: [ - C30004 ^2 + C30005 ^2 + C30006 ^2 ] <= 0
to get
qc0: [ - (0.5 +0.5 C0) ^2 + (32.46536616149585 C9003) ^2 + (0.5 - 0.5 C0) ^2 ] <= 0
which equals
qc0: - C0 + [ 1054 C9003 ^2 ] <= 0
These substitutions have to be made by hand. Thus I hope, that you can adjust your code in some way to achieve this reformulation easier.
0 -
The substitutions you mentioned can indeed be made on paper. These simplifications can be done because while implementing the cone constraints, some variables are repeated on both LHS and RHS. If I implement the cone constraint in its most simplified form, it involves multiplication of variables and this makes the individual terms non-linear. Hence, I reformulated the constraints so that they can be defined in the form of a standard second order cone.
From my understanding the constraint corresponding to the above lines of code is the last cone constraint in the PDF containing the equations I shared previously, i.e.
I did clean the code and improve the scaling. The LP file after that is on https://www.filemail.com/d/gmtqihzqiyxtpkm
Let me know if you observe any source of numerical issues. The maximum matrix coefficient range is 10^5 again. For reference, this is the output I have as of now:
+ Solver chosen : GUROBI-GUROBI
+ Processing objective function
+ Processing constraintsWarning: Some output might be missing due to a network interruption. To get the missing output, rerun the script.Set parameter FeasibilityTol to value 0.01
Set parameter OptimalityTol to value 0.01
Set parameter Method to value 2
Set parameter ScaleFlag to value 3
Set parameter BarConvTol to value 1
Set parameter BarHomogeneous to value 1
Set parameter CrossoverBasis to value 1
Set parameter NodefileDir to value ""
Set parameter NumericFocus to value 3
Set parameter Presolve to value 0
Set parameter PreSOS2BigM to value 0
Set parameter TuneTrials to value 3
Academic license - for non-commercial use only - expires 2022-06-18
Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 120011 rows, 114007 columns and 276013 nonzeros
Model fingerprint: 0x24a59e3a
Model has 18001 quadratic constraints
Coefficient statistics:
Matrix range [1e-03, 1e+02]
QMatrix range [1e+00, 1e+00]
Objective range [1e+00, 1e+00]
Bounds range [0e+00, 0e+00]
RHS range [8e-02, 3e+02]
Presolve time: 0.39s
Presolved: 120011 rows, 114007 columns, 276013 nonzeros
Presolved model has 18001 second-order cone constraints
Ordering time: 0.05s
Barrier statistics:
Free vars : 29995
AA' NZ : 3.660e+05
Factor NZ : 1.275e+06 (roughly 100 MB of memory)
Factor Ops : 1.479e+07 (less than 1 second per iteration)
Threads : 2
Objective Residual
Iter Primal Dual Primal Dual Compl Time
0 3.05989339e+02 -6.12063990e+05 4.69e+06 1.00e+03 9.39e+05 1s
1 4.30546591e+02 -4.05268237e+07 2.04e+06 6.63e+02 4.64e+05 3s
2 3.22184120e+02 -6.54174621e+07 4.50e+05 4.33e+02 1.71e+05 5s
3 2.97646221e+02 -9.08956509e+07 1.08e+05 1.12e+02 4.49e+04 8s
4 3.01067579e+02 -9.25209172e+06 3.07e+04 1.15e+01 2.81e+03 10s
5 3.04262863e+02 -2.16428337e+06 6.05e+03 2.61e+00 2.42e+02 11s
6 3.05110721e+02 -1.08123937e+06 1.45e+03 1.24e+00 6.49e+01 13s
7 3.05312817e+02 -4.09910578e+05 6.51e+02 4.57e-01 2.23e+01 15s
8 3.05350660e+02 -2.45910503e+05 5.03e+02 2.73e-01 1.34e+01 16s
9 3.05422626e+02 -2.05711025e+05 2.24e+02 2.28e-01 9.97e+00 18s
10 3.05429473e+02 -6.38419169e+04 1.36e+02 7.04e-02 3.20e+00 21s
11 3.05407242e+02 -3.72094206e+04 7.83e+01 4.11e-02 1.80e+00 23s
12 3.05350668e+02 -1.61275192e+04 5.04e+01 1.79e-02 8.08e-01 25s
13 3.05076437e+02 -1.02936272e+04 2.76e+01 1.14e-02 4.84e-01 27s
14 3.03610767e+02 -4.86884922e+03 1.14e+01 5.50e-03 2.12e-01 30s
15 3.01358629e+02 -1.24396262e+03 5.17e+00 1.62e-03 6.14e-02 33s
16 2.94026766e+02 -4.23471163e+02 1.39e+00 7.50e-04 2.58e-02 36s
17 2.86764753e+02 4.39243765e+00 8.08e-01 2.92e-04 1.04e-02 38s
18 2.69998359e+02 8.71365234e+01 5.38e-01 1.89e-04 7.18e-03 41s
19 2.40872670e+02 9.34706843e+01 4.41e-01 1.56e-04 6.66e-03 44s
20 1.81129080e+02 6.56029189e+01 3.44e-01 1.60e-04 7.29e-03 47s
21 1.12138041e+02 2.12475558e+01 3.33e-01 4.63e-04 8.61e-03 49s
22 9.10954179e+01 -1.49604268e+01 3.25e-01 3.55e-04 9.25e-03 50s
23 6.93762424e+01 -7.32996259e+01 3.15e-01 1.29e-03 9.92e-03 51s
24 4.93827849e+01 -1.09762021e+02 3.04e-01 3.25e-03 1.02e-02 54s
25 1.61263763e+00 -1.50688353e+02 2.85e-01 5.70e-03 1.08e-02 56s
26 -3.22854364e+01 -2.36473183e+02 2.74e-01 7.60e-03 1.21e-02 58s
27 -9.82573987e+01 -2.66430462e+02 2.51e-01 7.98e-03 1.28e-02 60s
28 -1.32649328e+02 -2.99170856e+02 2.43e-01 1.08e-02 1.33e-02 61s
29 -1.96800550e+02 -3.68987149e+02 2.28e-01 1.48e-02 1.45e-02 63s
30 -2.81997244e+02 -3.85271537e+02 2.14e-01 1.52e-02 1.55e-02 65s
31 -3.01538281e+02 -4.01694037e+02 2.14e-01 1.12e-02 1.59e-02 67s
32 -3.21125305e+02 -4.50860292e+02 2.08e-01 1.75e-02 1.64e-02 68s
33 -3.53660110e+02 -5.21310511e+02 2.04e-01 1.28e-02 1.74e-02 70s
34 -4.00877643e+02 -5.77749305e+02 1.97e-01 1.93e-02 1.83e-02 72s
35 -4.04558679e+02 -6.12549028e+02 1.96e-01 2.26e-02 1.86e-02 73s
36 -4.40930086e+02 -6.33383977e+02 1.90e-01 2.59e-02 1.90e-02 75s
37 -4.42663674e+02 -6.34496808e+02 1.88e-01 2.72e-02 1.88e-02 77s
38 -4.88961771e+02 -6.74144373e+02 1.83e-01 3.90e-02 1.95e-02 78s
39 -5.11303870e+02 -7.24564406e+02 1.80e-01 1.86e-02 2.03e-02 80s
40 -5.27345018e+02 -7.44977067e+02 1.76e-01 2.06e-02 2.05e-02 82s
41 -5.36353975e+02 -7.54015360e+02 1.71e-01 3.46e-02 2.02e-02 83s
42 -5.52270952e+02 -7.81983166e+02 1.68e-01 2.11e-02 2.04e-02 85s
43 -5.55750734e+02 -7.85276503e+02 1.66e-01 2.45e-02 2.03e-02 87s
44 -5.58528218e+02 -7.89303011e+02 1.65e-01 3.04e-02 2.02e-02 89s
45 -5.62566004e+02 -7.97749520e+02 1.64e-01 3.09e-02 2.02e-02 91s
46 -5.73634404e+02 -8.07535323e+02 1.61e-01 3.18e-02 2.03e-02 92s
47 -5.75300106e+02 -8.12332713e+02 1.61e-01 3.45e-02 2.03e-02 94s
48 -5.75605953e+02 -8.09851426e+02 1.60e-01 3.79e-02 2.02e-02 96s
49 -5.86065132e+02 -8.21579937e+02 1.59e-01 4.40e-02 2.03e-02 98s
50 -5.98101069e+02 -8.50965956e+02 1.57e-01 5.05e-02 2.07e-02 100s
51 -6.11145771e+02 -8.62918687e+02 1.56e-01 6.91e-02 2.09e-02 101s
52 -6.16248269e+02 -8.65317121e+02 1.54e-01 7.65e-02 2.08e-02 103s
53 -6.53533224e+02 -8.74810963e+02 1.49e-01 9.04e-02 2.10e-02 105s
54 -6.70231230e+02 -8.87906469e+02 1.46e-01 9.11e-02 2.11e-02 106s
55 -6.73254791e+02 -8.94683760e+02 1.45e-01 1.14e-01 2.11e-02 109s
56 -6.80399176e+02 -9.06256913e+02 1.44e-01 1.40e-01 2.12e-02 110s
57 -7.05816898e+02 -9.20461155e+02 1.42e-01 1.75e-01 2.14e-02 111s
58 -7.13859572e+02 -9.31281492e+02 1.41e-01 2.00e-01 2.16e-02 113s
59 -7.18082157e+02 -9.33976060e+02 1.41e-01 2.04e-01 2.16e-02 115s
60 -7.50275489e+02 -9.36204226e+02 1.38e-01 2.03e-01 2.18e-02 116s
61 -7.52636909e+02 -9.49673281e+02 1.38e-01 2.36e-01 2.19e-02 118s
62 -7.56023385e+02 -9.57730534e+02 1.38e-01 2.70e-01 2.19e-02 120s
63 -7.59913791e+02 -9.62141177e+02 1.37e-01 2.84e-01 2.20e-02 122s
64 -7.66713095e+02 -9.66384333e+02 1.37e-01 2.97e-01 2.21e-02 123s
65 -7.81501300e+02 -9.67190696e+02 1.34e-01 3.00e-01 2.21e-02 125s
66 -7.83069281e+02 -9.67888299e+02 1.34e-01 3.17e-01 2.20e-02 128s
67 -7.86902098e+02 -9.88333744e+02 1.34e-01 3.37e-01 2.22e-02 129s
68 -7.90377048e+02 -9.92025956e+02 1.33e-01 3.46e-01 2.23e-02 131s
69 -7.91531588e+02 -9.94989289e+02 1.33e-01 3.60e-01 2.23e-02 133s
70 -7.96655965e+02 -1.00450825e+03 1.32e-01 3.86e-01 2.24e-02 135s
Barrier performed 70 iterations in 134.84 seconds (13.37 work units)
Numerical trouble encounteredsol = struct with fields:
yalmipversion: '20210331'
matlabversion: '9.11.0.1769968 (R2021b)'
yalmiptime: 2.104350436266666e+02
solvertime: 3.082109563733333e+02
info: 'Numerical problems (GUROBI-GUROBI)'
problem: 40 -
Thank you for the files. We are still investigating. This might be a bug on our side but I cannot tell 100% right now and need to investigate further.
If I implement the cone constraint in its most simplified form, it involves multiplication of variables and this makes the individual terms non-linear.
Could you please try the nonconvex formulation? It might be possible that the spatial branch-and-bound algorithm is numerically more stable for this particular model. This is not something I would recommend in general but it may be helpful here.
0 -
Let me know the results of your investigation.
I did try the non-convex formulation prior to convex formulation. Unfortunately, the same numerical issues persisted.
0 -
Hi Tushar,
The numerical trouble come from the fact that very large values for barrier candidates are used. This is an indicator that your model is most likely unbounded (or infeasible is some cases).
I can't tell how to fix your model to be bounded. What definitely helps is to provide finite bounds for all variables, e.g., setting bounds of \(-10^5 \leq x \leq 10^5\) for all variables seemed to work a bit better.
Best regards,
Jaromił0 -
Hi Jeromil
Bounding the variables indeed solved the problem. The results are physically wrong, but that's the issue with problem formulation that I have to solve.
Could you elaborate on what barrier candidates are? Are they the same as optimisation variables? And how did you find them to be the source of the numerical issues? I simply checked the coefficients of the constraint variables in the LP file to improve the scaling. But if there's a means to conclude that the problem is unbounded from the LP file itself, I'm interested to learn that. From prior knowledge of the model, I know what values each variable would likely have. Since all of them are unique, I assumed the solution would converge to those values regardless of whether it's bounded or not.
Kind regards
Tushar Goyal0 -
Hi Tushar,
Could you elaborate on what barrier candidates are? Are they the same as optimisation variables? And how did you find them to be the source of the numerical issues?
A barrier candidate is a point in the solution space used to compute a step in the Barrier algorithm. These candidates are not available for the user. Thus, I had to look deeper and debug the code to get this information. The only indicator about this issue was the "Numerical trouble encountered" message but this message can have many sources.
I simply checked the coefficients of the constraint variables in the LP file to improve the scaling.
This is the correct way to address numerical issues as a user.
But if there's a means to conclude that the problem is unbounded from the LP file itself, I'm interested to learn that.
No, unfortunately not. A possible (and very difficult way) would be if you have enough knowledge about your model and are able to find some wrong constraint or a missing bound. This of course gets a lot harder with increased model size and complexity.
From prior knowledge of the model, I know what values each variable would likely have. Since all of them are unique, I assumed the solution would converge to those values regardless of whether it's bounded or not.
I think that this is very valuable knowledge and you should definitely use it to your advantage when modeling your problem. It is also very hard to think of all possibilities that may occur without providing all knowledge to the solver/model and if those lead to numerical difficulties, then spotting the source of the issue gets really difficult.
Bounding the variables indeed solved the problem. The results are physically wrong, but that's the issue with problem formulation that I have to solve.
Happy to hear that and I hope that further reformulation will improve the algorithm's behavior.
Best regards,
Jaromił0
Please sign in to leave a comment.
Comments
15 comments