python - 使用 CVXPY 求解拟凸问题

Question

我正在尝试使用 CVXPY 包(使用 ECOS 求解器)解决 Python 中的优化问题。该问题涉及最小化受多个约束的线性变量。问题本身是拟凸的，所以我实现了二分算法来找到最优解。

问题在于，当二分越来越接近最优值时，求解器总是失败。这就像当解决方案接近最优时，求解器不知道该怎么做。我相信这个问题与容差问题有关，但我不确定要更改哪个容差变量。

我尝试在 Python 3.6 中编写这个问题。我也尝试过使用不同的求解器，但求解器在接近最优时总是失败。代码如下所示。请注意，“gamma”是我在这里最小化的变量。

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)

这是输出。您可以看到解收敛到接近 1.059...但随后解算器遇到错误。

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]

这是最后一次二分迭代的结果 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 1

问题可能出在这里:

cp.multiply(cp.inv_pos(gamma),cp.abs(cp.multiply(Wd,PSI - cp.multiply(G,To)))) - cp.real(PSI)

cp.inv_pos(gamma) 是凸非负函数，cp.abs(...) 也是如此。

当参数非负时，标量积是拟凹的，就像这里的情况一样。但为了满足拟凹性，两个输入都必须是凹的。由于输入是凸的，因此曲率未知。因此，问题本身是非凸的。 (参见here。)