NUMERICAL PROBLEMS / Slacks/multipliers leaving the cone
Closed this issue · 3 comments
dset0x commented
Hello,
I have the following program:
#!/usr/bin/env python
import cvxpy as cp
from cvxpy.settings import OPTIMAL, OPTIMAL_INACCURATE
from math import sqrt
from typing import *
Vbase = 231
Sbase = 1e6
Zbase = Vbase**2 / Sbase
buses: List[dict] = [
{'id':1, 'pc':0, 'pg':0, 'qc':0, 'sg':0},
{'id':2, 'pc':5e03, 'pg':3e03, 'qc':5e03, 'sg':10e03},
{'id':3, 'pc':4e03, 'pg':3e03, 'qc':6.113e03, 'sg':10e03},
{'id':4, 'pc':10e03, 'pg':3e03, 'qc':0e03, 'sg':10e03},
]
lines = [
{'r':0.09323, 'x':0.05585},
{'r':0.19495, 'x':0.0812},
{'r':0.11738, 'x':0.02925},
]
V0=231.0
epsilon=0.05
##############################################################################
V0 = V0/Vbase
r = [line['r']/Zbase for line in lines]
x = [line['x']/Zbase for line in lines]
Q = cp.Variable((len(buses)))
P = cp.Variable((len(buses)))
U = cp.Variable((len(buses)))
pc = [bus['pc']/Sbase for bus in buses]
qc = [bus['qc']/Sbase for bus in buses]
qg = [cp.Parameter(value=0), cp.Variable(), cp.Variable(), cp.Variable()]
pg = [bus['pg']/Sbase for bus in buses]
sg = [bus['sg']/Sbase for bus in buses]
s_h = [sqrt(sg[j]**2 - pg[j]**2) for j in range(len(buses))]
constraints = []
for j in range(len(buses)-1):
jp1 = j+1
constraints.extend((
P[jp1] == (P[j] - pc[jp1] + pg[jp1]),
Q[jp1] == (Q[j] - qc[jp1] + qg[jp1]),
U[jp1] == (U[j] - 2 * (r[j] * P[j] + x[j] * Q[j])),
))
constraints.extend((
P[len(buses)-1] == 0,
Q[len(buses)-1] == 0,
U[0] == 0,
))
for j in range(1, len(buses)):
constraints.extend((
V0**2 * (epsilon**2 - 2*epsilon) <= U[j],
U[j] <= V0**2 * (epsilon**2 + 2*epsilon),
cp.abs(qg[j]) <= s_h[j],
))
prob = cp.Problem(
cp.Minimize(
cp.sum(
cp.multiply(
r,
(P[:-1]**2 + Q[:-1]**2) / V0**2
)
)
),
constraints
)
z=cp.sum(
cp.multiply(
r,
(P[:-1]**2 + Q[:-1]**2) / V0**2
)
)
prob.solve(verbose=True, solver=cp.ECOS)
if not prob.status in (OPTIMAL, OPTIMAL_INACCURATE):
raise Exception(prob.status)
assert 6600 < qg[1].value * Sbase < 6900
assert 7600 < qg[2].value * Sbase < 7800
assert 2200 < qg[3].value * Sbase < 2500running under the following environment:
pyenv install 3.6.7
pyenv virtualenv 3.6.7 foo
pyenv activate foo
pip --no-cache-dir install cvxpyWhich fails with
ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
It pcost dcost gap pres dres k/t mu step sigma IR | BT
0 +3.888e-04 -3.519e+03 +2e-201 inf 0e+00 7e-01 8e-203 --- ---Slacks/multipliers leaving the cone, recovering best iterate (0) and stopping.
NUMERICAL PROBLEMS (reached feastol=inf, reltol=-nan, abstol=1.7e-201).Traceback (most recent call last):
File "/home/user/.pyenv/versions/3.6.7/envs/foo/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 635, in unpack_results
self.unpack(solution)
File "/home/user/.pyenv/versions/3.6.7/envs/foo/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 610, in unpack
raise ValueError("Cannot unpack invalid solution.")
ValueError: Cannot unpack invalid solution.
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "./bcos.py", line 65, in <module>
prob.solve(verbose=True, solver=cp.ECOS)
File "/home/user/.pyenv/versions/3.6.7/envs/foo/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 271, in solve
return solve_func(self, *args, **kwargs)
File "/home/user/.pyenv/versions/3.6.7/envs/foo/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 506, in _solve
self.unpack_results(solution, full_chain, inverse_data)
File "/home/user/.pyenv/versions/3.6.7/envs/foo/lib/python3.6/site-packages/cvxpy/problems/problem.py", line 639, in unpack_results
"Try another solver or solve with verbose=True for more "
cvxpy.error.SolverError: Solver 'ECOS' failed. Try another solver or solve with verbose=True for more information. Try recentering the problem data around 0 and rescaling to reduce the dynamic range.
In a different environment (with python 3.7.x) I was getting different results every time I executed it (either failure like above or very bad results).
I have a colleague for whom it seems to be working, but I have not dug into the differences between the systems yet (I know the python, ecos and cvxpy versions are the same).
It works okay with OSQP. Are you perhaps able to reproduce this?
dset0x commented
I just noticed the superfluous block
z=cp.sum(
cp.multiply(
r,
(P[:-1]**2 + Q[:-1]**2) / V0**2
)
)Removing it might get you out of the exception and into bad results. Or putting it back in might. Or it might not. Sometimes adding an import time; time.sleep(2) before solve() might trigger this, or not.
dset0x commented
Setting the non-determinism of yesterday aside, I've created the same setup on two other machines and, indeed, ECOS runs successfully there. So far we have:
- Machine A: Gentoo x86_64 on i7-4500U - Fails (With and without latest ucode update)
ipdb> pl.architecture(), pl.machine(), pl.platform(), pl.python_build(), pl.python_compiler(), pl.python_implementation(), pl.python_version(), pl.release(), pl.system(), pl.version(),
(('64bit', 'ELF'), 'x86_64', 'Linux-4.19.27-gentoo-r1-x86_64-Intel-R-_Core-TM-_i7-4500U_CPU_@_1.80GHz-with-gentoo-2.6', ('default', 'Apr 10 2019 09:05:50'), 'GCC 7.3.0', 'CPython', '3.6.7', '4.19.27-gentoo-r1', 'Linux', '#1 SMP Sat Apr 6 11:30:35 CEST 2019')
- Machine B: Gentoo x86_64 on Celeron E3300 - Works
ipdb> pl.architecture(), pl.machine(), pl.platform(), pl.python_build(), pl.python_compiler(), pl.python_implementation(), pl.python_version(), pl.release(), pl.system(), pl.version(),
(('64bit', ''), 'x86_64', 'Linux-4.14.83-gentoo-x86_64-Intel-R-_Celeron-R-_CPU_E3300_@_2.50GHz-with-gentoo-2.6', ('default', 'Apr 10 2019 09:03:31'), 'GCC 7.3.0', 'CPython', '3.6.7', '4.14.83-gentoo', 'Linux', '#1 SMP Fri Jan 4 16:51:26 EET 2019')
- Machine C: Devuan i386 on T2130 - Works
Stepping though cvxpy in the code above on machines A and B we land at /home/user/.pyenv/versions/foo/lib/python3.6/site-packages/cvxpy/reductions/solvers/conic_solvers/ecos_conif.py(148)solve_via_data(). Both result in the same data:
> print(f"C:\n{data[s.C]}\n\nG:\n{data[s.G]}\n\nH:\n{data[s.H]}\n\nCones:\n{cones}\n\nA:\n{data[s.A]}\n\nB:\n{data[s.B]}")
C:
[1.74715616 3.65341729 2.19973389 1.74715616 3.65341729 2.19973389
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. ]
G:
(15, 0) -1.0
(16, 0) -1.0
(18, 1) -1.0
(19, 1) -1.0
(21, 2) -1.0
(22, 2) -1.0
(24, 3) -1.0
(25, 3) -1.0
(27, 4) -1.0
(28, 4) -1.0
(30, 5) -1.0
(31, 5) -1.0
(17, 6) -2.0
(20, 7) -2.0
(23, 8) -2.0
(26, 10) -2.0
(29, 11) -2.0
(32, 12) -2.0
(2, 14) 1.0
(3, 14) -1.0
(0, 16) -1.0
(1, 16) 1.0
(5, 17) -1.0
(6, 17) 1.0
(10, 18) -1.0
(11, 18) 1.0
(7, 19) 1.0
(8, 19) -1.0
(12, 20) 1.0
(13, 20) -1.0
(2, 21) -1.0
(3, 21) -1.0
(4, 21) 1.0
(7, 22) -1.0
(8, 22) -1.0
(9, 22) 1.0
(12, 23) -1.0
(13, 23) -1.0
(14, 23) 1.0
H:
[ 0.0975 0.1025 -0. -0. 0.00953939 0.0975
0.1025 -0. -0. 0.00953939 0.0975 0.1025
-0. -0. 0.00953939 1. -1. 0.
1. -1. 0. 1. -1. 0.
1. -1. 0. 1. -1. 0.
1. -1. 0. ]
Cones:
{'l': 15, 'q': [3, 3, 3, 3, 3, 3], 'e': 0}
A:
(0, 6) -1.0
(2, 6) 3.4943123254811566
(0, 7) 1.0
(3, 7) -1.0
(5, 7) 7.306834579561853
(3, 8) 1.0
(6, 8) -1.0
(8, 8) 4.399467776091153
(6, 9) 1.0
(9, 9) 1.0
(1, 10) -1.0
(2, 10) 2.0932891062761194
(1, 11) 1.0
(4, 11) -1.0
(5, 11) 3.043421225239407
(4, 12) 1.0
(7, 12) -1.0
(8, 12) 1.096306291111486
(7, 13) 1.0
(10, 13) 1.0
(1, 14) -1.0
(2, 15) -1.0
(11, 15) 1.0
(2, 16) 1.0
(5, 16) -1.0
(5, 17) 1.0
(8, 17) -1.0
(8, 18) 1.0
(4, 19) -1.0
(7, 20) -1.0
B:
[-0.002 -0.005 -0. -0.001 -0.006113 -0. -0.007
-0. -0. -0. -0. -0. ]
Therefore numpy seems to be working the same on both machines, while ECOS does not:
Machine A:
ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
It pcost dcost gap pres dres k/t mu step sigma IR | BT
0 +9.702e+00 +2.298e-01 +0e+00 0e+00 0e+00 1e+00 0e+00 --- --- 0 0 - | - -
OPTIMAL (within feastol=0.0e+00, reltol=0.0e+00, abstol=0.0e+00).
Runtime: 0.000204 seconds.
--------------------------------------------------------------------------------------------------------------------------
ipdb> pprint(solution)
{'info': {'dcost': 0.22982938543193038,
'dinf': 0.0,
'dinfres': nan,
'dres': 0.0,
'exitFlag': 0,
'gap': 0.0,
'infostring': 'Optimal solution found',
'iter': 0,
'mi_iter': -1,
'numerr': 0,
'pcost': 9.701526081563065,
'pinf': 0.0,
'pinfres': nan,
'pres': 0.0,
'r0': 1e-08,
'relgap': 0.0,
'timing': {'runtime': 0.000204425,
'tsetup': 0.000128833,
'tsolve': 7.5592e-05}},
's': array([inf, nan, nan, nan, inf, inf, inf, inf, inf, inf, inf, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, inf, inf, inf, inf,
inf, inf, inf, inf, inf, inf, inf]),
'x': array([1.26287353e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
1.27067061e+00, 1.29688362e+00, 3.60896340e-01, 3.81579800e-13,
2.00414011e-06, 9.99972088e-01, 5.65065001e-01, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 1.57401079e+59, 3.37510073e-04]),
'y': array([1.07850946e+00, 0.00000000e+00, 1.61806591e-03, 1.70398946e-03,
3.20797645e-03, 1.14155965e-03, 3.39128939e-03, 9.77669014e-04,
1.52943527e-03, 3.08576782e-03, 3.54405017e-03, 3.48294586e-03]),
'z': array([6.74626390e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
6.74626388e-01, 6.74626388e-01, 6.74626388e-01, 9.93339630e-01,
2.79092443e-05, 9.99972088e-01, 6.74626385e-01, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.33421652e-01,
2.33421652e-01, 2.22516953e-01, 4.50487776e-13, 7.50812960e-14,
1.99429999e-01, 2.95615437e-01, 4.71039344e-13, 1.04504913e-13,
2.06479420e-01, 2.06479424e-01, 0.00000000e+00, 1.96509831e-01,
5.63601168e-13])}
Machine B:
ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS
It pcost dcost gap pres dres k/t mu step sigma IR | BT
0 +1.204e-16 -2.519e+00 +5e+01 7e-01 8e-02 1e+00 2e+00 --- --- 1 1 - | - -
1 -8.302e-01 -9.302e-01 +3e+00 7e-02 3e-03 6e-02 1e-01 0.9603 1e-02 1 1 1 | 0 0
2 -5.039e-02 -5.689e-02 +2e-01 4e-03 2e-04 3e-03 9e-03 0.9301 3e-03 2 2 2 | 0 0
3 -4.068e-03 -4.747e-03 +2e-02 3e-04 1e-05 2e-04 9e-04 0.9315 4e-02 1 2 2 | 0 0
4 -4.761e-04 -6.526e-04 +5e-03 9e-05 4e-06 5e-05 2e-04 0.7517 3e-02 1 1 2 | 0 0
5 -1.397e-04 -2.416e-04 +3e-03 6e-05 2e-06 4e-05 2e-04 0.6099 4e-01 1 1 1 | 0 0
6 +5.234e-04 +5.076e-04 +5e-04 9e-06 4e-07 6e-06 2e-05 0.8529 1e-02 1 1 1 | 0 0
7 +6.162e-04 +6.129e-04 +1e-04 2e-06 9e-08 2e-06 6e-06 0.9890 2e-01 1 1 1 | 0 0
8 +6.429e-04 +6.428e-04 +4e-06 6e-08 3e-09 6e-08 2e-07 0.9701 6e-04 2 1 1 | 0 0
9 +6.438e-04 +6.438e-04 +8e-08 1e-09 6e-11 1e-09 4e-09 0.9777 1e-04 1 1 1 | 0 0
10 +6.438e-04 +6.438e-04 +4e-09 7e-11 3e-12 6e-11 2e-10 0.9502 1e-04 1 1 1 | 0 0
OPTIMAL (within feastol=6.9e-11, reltol=6.4e-06, abstol=4.1e-09).
Runtime: 0.005448 seconds.
--------------------------------------------------------------------------------------------------------------------------
ipdb> pprint(solution)
{'info': {'dcost': 0.0006438011394196911,
'dinf': 0.0,
'dinfres': nan,
'dres': 2.9100825007592447e-12,
'exitFlag': 0,
'gap': 4.0952666110284685e-09,
'infostring': 'Optimal solution found',
'iter': 10,
'mi_iter': -1,
'numerr': 0,
'pcost': 0.0006438012511756154,
'pinf': 0.0,
'pinfres': 0.816875053174683,
'pres': 6.854394881079808e-11,
'r0': 1e-08,
'relgap': 6.361073878682253e-06,
'timing': {'runtime': 0.00544805,
'tsetup': 0.000148619,
'tsolve': 0.005299431}},
's': array([ 7.45321338e-02, 1.25467866e-01, 1.63584218e-03, 1.51193082e-02,
1.16181695e-03, 2.81873888e-02, 1.71812611e-01, 8.94485316e-04,
1.63191683e-02, 9.32565301e-04, 6.02644231e-09, 1.99999994e-01,
6.60892541e-03, 1.13683465e-02, 5.50756144e-04, 1.00010000e+00,
-9.99900000e-01, 2.00000000e-02, 1.00006400e+00, -9.99936000e-01,
1.60000000e-02, 1.00004900e+00, -9.99951000e-01, 1.40000000e-02,
1.00003273e+00, -9.99967273e-01, -1.14415701e-02, 1.00001583e+00,
-9.99984167e-01, -7.95810412e-03, 1.00000566e+00, -9.99994337e-01,
-4.75942112e-03]),
'x': array([ 9.99999397e-05, 6.39999358e-05, 4.89999380e-05, 3.27272977e-05,
1.58327859e-05, 5.66294417e-06, 9.99999999e-03, 7.99999999e-03,
6.99999999e-03, -4.24740487e-12, -5.72078504e-03, -3.97905206e-03,
-2.37971056e-03, 8.95019745e-14, 6.74173298e-03, 2.70852792e-13,
-2.29678662e-02, -6.93126113e-02, -9.74999940e-02, 7.71234150e-03,
2.37971056e-03, 8.37757512e-03, 8.60682676e-03, 8.98863592e-03]),
'y': array([ 6.83124438e-02, 6.30362358e-08, 9.54961239e-03, 1.96544409e-01,
5.95711919e-08, 9.54961207e-03, 2.69353874e-01, 1.05583884e-08,
9.54961021e-03, -2.69353874e-01, -1.05583884e-08, 9.54961239e-03]),
'z': array([ 8.25643751e-10, 5.04223239e-10, 6.70121778e-08, 3.97639434e-09,
7.09389501e-08, 2.23767159e-09, 3.79234957e-10, 6.37037705e-08,
4.13238745e-09, 6.77868078e-08, 9.54961053e-03, 3.20778864e-10,
1.72446811e-08, 6.68629041e-09, 2.38815441e-08, 8.73665439e-01,
8.73490724e-01, -1.74715578e-02, 1.82682555e+00, 1.82659174e+00,
-2.92272649e-02, 1.09992084e+00, 1.09981305e+00, -1.53981310e-02,
8.73606671e-01, 8.73549492e-01, 9.99501828e-03, 1.82673755e+00,
1.82667974e+00, 1.45317478e-02, 1.09987317e+00, 1.09986072e+00,
5.23467338e-03])}
dset0x commented
Not a bug in ecos after all -- see issue referenced above.
stderred was overriding the init function of ecos, meaning w was never initialized and its contents were left at the mercy of malloc. Explains the non-deterministic behavior.