Нахождение решения квазивыпуклой задачи с использованием CVXPY - PullRequest
Новая
       46

Нахождение решения квазивыпуклой задачи с использованием CVXPY

1 голос
/ 05 июня 2019

Я пытаюсь решить проблему оптимизации в Python с помощью пакета CVXPY (с решателем ECOS). Проблема заключается в минимизации линейной переменной с учетом множества ограничений. Сама проблема квазивыпуклая, поэтому я реализую алгоритм деления пополам, чтобы найти оптимальное решение.

Проблема в том, что решатель всегда терпит неудачу, когда деление пополам приближается к оптимальному значению. Как будто решатель не знает, что делать, когда решение приближается к оптимальному. Я считаю, что эта проблема связана с вопросом толерантности, но я не уверен, какую переменную допуска изменить.

Я пытался кодировать эту проблему в Python 3.6. Я также пытался использовать разные решатели, но решатель всегда терпит неудачу при приближении к оптимальному. Код показан ниже. Обратите внимание, что «гамма» - это переменная, которую я здесь минимизирую. 

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.```

1 Ответ

0 голосов
/ 11 июля 2019

Проблема может быть здесь: 

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(...).

Скалярные произведения являются квазивогнутыми, когда их аргументы неотрицательны, как в данном случае. Но чтобы удовлетворить квазивогнутость, оба входа должны быть вогнутыми. Поскольку входы выпуклые, кривизна неизвестна. Поэтому сама проблема невыпуклая. (См. здесь .)

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...