Finding solution to quasi-convex problem using CVXPY
An answer to this question on Stack Overflow.
Question
I am trying to solve an optimization problem in Python using the CVXPY package (with the ECOS solver). The problem involves minimizing a linear variable subject to multiple constraints. The problem itself is quasi-convex, so I implement a bisection algorithm to find the optimal solution.
The problem is that the solver always fails when the bisection is getting closer and closer to the optimal value. It's like the solver doesn't know what to do when the solution is approaching the optimal. I believe this problem is related to a tolerance issue, but I'm not sure what tolerance variable to change.
I have tried coding this problem in Python 3.6. I have also tried using different solvers, but the solver always fails when approaching the optimal. The code is shown below. Note that 'gamma' is the variable I am minimizing here.
import cvxpy as cp
import numpy as np
Ts = 500e-6;
w = np.logspace(-1, np.log10(np.pi/Ts), 400);
R = 0.1; L = 4e-2; a = R/L;
G = 0.0125/(np.exp(1j*w*Ts) - np.exp(-a*Ts))
bw = 200; zeta = 0.8;
wn = (2*np.pi*bw)/np.sqrt(1 - 2*np.power(zeta,2) + np.sqrt(2 - 4*np.power(zeta,2) + 4*np.power(zeta,4)))
Td = np.power(wn,2)/(-np.power(w,2) + 2*zeta*wn*1j*w + np.power(wn,2)); Wd = 1/(1-Td);
n = 5
rho_r = cp.Variable(n); rho_s = cp.Variable(n-1); rho_t = cp.Variable(n);
Ro = 0;
for i in range(n):
Ri = rho_r[i]*np.exp(-i*1j*w*Ts);
Ro = Ro + Ri;
So = 0
for i in range(n-1):
Si = rho_s[i]*np.exp(-(i+1)*1j*w*Ts);
So = So + Si;
So_no_int = 1+So
So = cp.multiply((1-np.exp(-1j*w*Ts)),(1+So));
To = 0;
for i in range(n):
Ti = rho_t[i]*np.exp(-i*1j*w*Ts);
To = To + Ti;
g_max = 2; g_min = 1e-8; g_tol = 1e-5; gamma = (g_max + g_min)/2;
mm = 0.5;
while g_max - g_min > g_tol:
PSI = So + cp.multiply(G,Ro)
F1 = cp.multiply(cp.inv_pos(gamma),cp.abs(cp.multiply(Wd,PSI - cp.multiply(G,To)))) - cp.real(PSI)
F2 = cp.abs(cp.multiply(mm,So)) - cp.real(PSI)
constraints = [F1 <= -1e-6, F2 <= -1e-6, cp.real(So_no_int) >= 5e-3]
prob = cp.Problem(cp.Minimize(0),constraints)
prob.solve(solver = cp.ECOS,feastol_inacc = 1e-7)
stat = prob.status
if stat == "optimal":
print("Feasible (gamma = ", gamma ,")")
gamma_opt = gamma;
rho_r_OPT = rho_r.value;
rho_s_OPT = rho_s.value;
rho_t_OPT = rho_t.value;
g_max = gamma;
gamma = np.average([gamma,g_min]);
GAIN = np.sum(rho_t_OPT)/np.sum(rho_r_OPT);
else:
print("Infeasible (gamma = ", gamma ,")")
g_min = gamma;
gamma = np.average([gamma,g_max]);
print("\nThe optimal solution gamma = %f" % gamma_opt)
Here is the output. You can see that the solution is converging to something close to 1.059....but then the solver runs into an error.
Infeasible (gamma = 1.000000005 )
Feasible (gamma = 1.5000000025 )
Feasible (gamma = 1.2500000037499999 )
Feasible (gamma = 1.125000004375 )
Feasible (gamma = 1.0625000046875 )
Infeasible (gamma = 1.0312500048437498 )
Infeasible (gamma = 1.0468750047656248 )
Infeasible (gamma = 1.0546875047265623 )
Infeasible (gamma = 1.0585937547070312 )
Feasible (gamma = 1.0605468796972657 )
Feasible (gamma = 1.0595703172021484 )
Infeasible (gamma = 1.0590820359545898 )
Traceback (most recent call last):
File "X:\Google Drive\...
Stuff\...\Python_Controller_optimization\main_no_LMIs.py", line 60, in <module>
prob.solve(solver = cp.ECOS,feastol_inacc = 1e-7)
File "C:\Users\...\AppData\Local\Programs\Python\Python37-32\lib\site-packages\cvxpy\problems\problem.py", line 289, in solve
return solve_func(self, *args, **kwargs)
File "C:\Users\...\AppData\Local\Programs\Python\Python37-32\lib\site-packages\cvxpy\problems\problem.py", line 574, in _solve
self.unpack_results(solution, full_chain, inverse_data)
File "C:\Users\...\AppData\Local\Programs\Python\Python37-32\lib\site-packages\cvxpy\problems\problem.py", line 717, 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.
[Finished in 4.8s]
Here is the result of the last bisection iteration with verbose = True:
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 +0.000e+00 -6.380e+02 +3e+03 1e-02 1e+01 1e+00 2e+00 --- --- 1 0 - | - -
1 +0.000e+00 -1.145e+02 +8e+02 4e-04 4e+00 1e+00 4e-01 0.8673 1e-01 1 1 1 | 0 0
2 +0.000e+00 -7.961e+01 +4e+02 3e-04 2e+00 6e-01 2e-01 0.7724 4e-01 2 1 2 | 0 0
3 +0.000e+00 -3.442e+01 +2e+02 1e-04 8e-01 2e-01 9e-02 0.6591 2e-01 2 1 2 | 0 0
4 +0.000e+00 -1.446e+01 +1e+02 6e-05 2e-01 5e-02 5e-02 0.7673 4e-01 2 2 2 | 0 0
5 +0.000e+00 -1.551e+00 +1e+01 6e-06 2e-02 4e-03 5e-03 0.8955 2e-02 2 1 1 | 0 0
6 +0.000e+00 -7.516e-01 +5e+00 3e-06 1e-02 1e-03 3e-03 0.7389 3e-01 2 2 2 | 0 0
7 +0.000e+00 -2.738e-01 +2e+00 1e-06 3e-03 5e-04 1e-03 0.7857 2e-01 3 3 3 | 0 0
8 +0.000e+00 -1.707e-01 +1e+00 7e-07 1e-03 3e-04 6e-04 0.7419 5e-01 3 3 3 | 0 0
9 +0.000e+00 -3.000e-02 +2e-01 1e-07 2e-04 4e-05 1e-04 0.8592 4e-02 2 3 3 | 0 0
10 +0.000e+00 -2.762e-02 +2e-01 1e-07 1e-04 4e-05 1e-04 0.1937 6e-01 3 3 4 | 0 0
11 +0.000e+00 -1.464e-02 +1e-01 6e-08 4e-05 2e-05 6e-05 0.6035 2e-01 3 4 4 | 0 0
12 +0.000e+00 -1.150e-02 +1e-01 5e-08 1e-05 2e-05 5e-05 0.5393 6e-01 3 4 4 | 0 0
13 +0.000e+00 -8.110e-03 +1e-01 3e-08 5e-06 1e-05 5e-05 0.4090 3e-01 4 5 5 | 0 0
14 +0.000e+00 -4.994e-03 +9e-02 2e-08 2e-06 8e-06 5e-05 0.5449 3e-01 5 6 6 | 0 0
15 +0.000e+00 -4.390e-03 +1e-01 2e-08 2e-06 8e-06 5e-05 0.2976 6e-01 5 6 6 | 0 0
16 +0.000e+00 -1.789e-03 +6e-02 8e-09 5e-07 4e-06 3e-05 0.6651 1e-01 4 5 5 | 0 0
17 +0.000e+00 -9.008e-04 +4e-02 5e-09 2e-07 2e-06 2e-05 0.5673 1e-01 5 6 6 | 0 0
18 +0.000e+00 -6.047e-05 +9e-03 3e-09 1e-08 4e-07 4e-06 0.9890 6e-02 5 5 6 | 0 0
19 +0.000e+00 -1.067e-05 +3e-03 3e-09 4e-09 3e-07 1e-06 0.9708 2e-01 5 6 6 | 0 0
20 +0.000e+00 -8.619e-07 +3e-04 3e-09 4e-09 3e-08 2e-07 0.9503 3e-02 6 6 6 | 0 0
21 +0.000e+00 -1.358e-07 +6e-05 3e-09 4e-09 6e-09 3e-08 0.8718 3e-02 6 5 6 | 0 0
22 +0.000e+00 -1.090e-07 +5e-05 3e-09 4e-09 7e-09 3e-08 0.4582 6e-01 5 5 5 | 0 0
23 +0.000e+00 -3.494e-08 +2e-05 3e-09 4e-09 4e-09 1e-08 0.8998 2e-01 7 5 5 | 0 0
24 +0.000e+00 -2.603e-08 +2e-05 3e-09 4e-09 3e-09 9e-09 0.5242 5e-01 7 5 5 | 0 0
25 +0.000e+00 -2.509e-08 +2e-05 3e-09 4e-09 3e-09 9e-09 0.1730 8e-01 7 5 5 | 0 0
26 +0.000e+00 -1.781e-08 +1e-05 3e-09 4e-09 3e-09 7e-09 0.4828 4e-01 5 5 5 | 0 0
27 +0.000e+00 -1.782e-08 +1e-05 3e-09 4e-09 3e-09 7e-09 0.0071 1e+00 6 5 5 | 0 0
28 +0.000e+00 -7.540e-09 +6e-06 3e-09 4e-09 2e-09 3e-09 0.9672 4e-01 5 5 6 | 0 0
29 +0.000e+00 -7.547e-09 +6e-06 3e-09 4e-09 2e-09 3e-09 0.0134 9e-01 5 5 5 | 0 0
30 +0.000e+00 -7.557e-09 +8e-06 1e-05 4e-09 2e-09 4e-09 0.1284 1e+00 0 4 6 | 0 0
Unreliable search direction detected, recovering best iterate (18) and stopping.
NUMERICAL PROBLEMS (reached feastol=1.1e-08, reltol=-nan(ind), abstol=8.9e-03).
Runtime: 0.182235 seconds.```
# Answer
The problem might be here:
cp.multiply(cp.inv_pos(gamma),cp.abs(cp.multiply(Wd,PSI - cp.multiply(G,To)))) - cp.real(PSI)
`cp.inv_pos(gamma)` is a convex, non-negative function and so is `cp.abs(...)`.
Scalar products are quasiconcave when their arguments are non-negative, as is the case here. But to satisfy quasiconcavity, both of the inputs must be concave. Since the inputs are convex, the curvature is not known. Therefore, the problem itself is non-convex. (See [here](https://www.cvxpy.org/tutorial/dqcp/index.html).)