Конвертировать полуопределенную программу из CVX в CVXPY

Question

Я хочу преобразовать следующий SDP - который только проверяет выполнимость ограничений - из CVX (MATLAB) в CVXPY (Python):

Ah = [1.0058, -0.0058; 1, 0];
Bh = [-1; 0];
Ch = [1.0058, -0.0058; -0.9829, 0.0056];
Dh = [-1; 1];

M = [0, 1;1, 0];
ni = size(M,1)/2;
n = size(Ah,1);
rho = 0.5;

cvx_begin sdp quiet
    variable P(n,n) semidefinite
    variable lambda(ni) nonnegative
    Mblk = M*kron(diag(lambda),eye(2));
    lambda(ni) == 1  % break homogeneity (many ways to do this...)
    [Ah Bh]'*P*[Ah Bh] - rho^2*blkdiag(P,0) + [Ch Dh]'*Mblk*[Ch Dh] <= 0
cvx_end


switch cvx_status
    case 'Solved'
        feas = 1;
    otherwise
        feas = 0;
end

Ниже мой код Python,

import cvxpy as cvx
import numpy as np
import scipy as sp


Ah = np.array([[1.0058, -0.0058], [1, 0]])
Bh = np.array([[-1], [0]])
Ch = np.array([[1.0058, -0.0058], [-0.9829, 0.0056]])
Dh = np.array([[-1], [1]])

M = np.array([[0, 1], [1, 0]])
ni, n = M.shape[0] / 2, Ah.shape[0]
rho = 0.5

P = cvx.Semidef(n)
lamda = cvx.Variable()

Mblk = np.dot(M, np.kron(cvx.diag(lamda), np.eye(2)))
ABh = np.concatenate((Ah, Bh), axis=1)
CDh = np.concatenate((Ch, Dh), axis=1)
constraints = [lamda[-1] == 1,
               np.dot(ABh.T, np.dot(P, ABh)) - rho**2*np.linalg.block_diag(P, 0) +
               np.dot(CDh.T, np.dot(Mblk, CDh)) << 0]

prob = cvx.Problem(cvx.Minimize(1), constraints)
feas = prob.status is cvx.OPTIMAL

Есть несколько ошибок при запуске программы. 1. Когда я печатаю Mblk, он показывает

Traceback (последний вызов был последним):

Файл "/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py", строка 2820, в run_code

Out [1]: exec code_obj в self.user_global_ns, self.user_ns

Файл "", строка 1, в

Mblk

Файл "/usr/lib/python2.7/dist-packages/IPython/core/displayhook.py", строка 247, в вызове

format_dict, md_dict = self.compute_format_data(результат)

Файл "/usr/lib/python2.7/dist-packages/IPython/core/displayhook.py", строка 157, в compute_format_data

вернуть self.shell.display_formatter.format(результат)

Файл "/usr/lib/python2.7/dist-packages/IPython/core/formatters.py", строка 152, в формате

data = formatter(obj)

Файл "/usr/lib/python2.7/dist-packages/IPython/core/formatters.py", строка 481, в вызове

printer.pretty(OBJ)

Файл "/usr/lib/python2.7/dist-packages/IPython/lib/pretty.py", строка 362, в довольно

вернуть _default_pprint(объект, себя, цикл)

Файл "/usr/lib/python2.7/dist-packages/IPython/lib/pretty.py", строка 482, в _default_pprint

p.text (магнезии (OBJ))

Файл "/usr/lib/python2.7/dist-packages/numpy/core/numeric.py", строка 1553, в array_repr

',', "array (")

Файл "/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py", строка 454, в array2string

разделитель, префикс, formatter=formatter)

Файл "/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py", строка 256, в _array2string

'int': IntegerFormat (данные),

Файл "/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py", строка 641, в init

max_str_len = max (len (str (Maximum.reduce(data))),

Файл "/usr/local/lib/python2.7/dist-packages/cvxpy/constraints/leq_constraint.py", строка 67, ненулевой

Возникновение исключения ("Невозможно оценить значение истинности ограничения.")

Исключение: Невозможно оценить истинное значение ограничения.

Когда я подойду к этой линии,

  constraints = [lamda[-1] == 1,
                   np.dot(ABh.T, np.dot(P, ABh)) - rho**2*np.linalg.block_diag(P, 0) +
                   np.dot(CDh.T, np.dot(Mblk, CDh)) << 0]

это показывает

Traceback (последний вызов был последним): файл

"... / sdp.py", строка 22, в

np.dot(ABh.T, np.dot(P, ABh)) - rho**2*np.linalg.block_diag(P, 0) + 

ValueError: установка элемента массива с последовательностью.

Как исправить эти проблемы?

Answer 1

Большая проблема с вашим кодом заключается в том, что вы не можете использовать функции NumPy для объектов CVXPY. Вам необходимо использовать эквивалентные функции CVXPY. Вот рабочая версия вашего кода:

import cvxpy as cvx
import numpy as np
import scipy as sp


Ah = np.array([[1.0058, -0.0058], [1, 0]])
Bh = np.array([[-1], [0]])
Ch = np.array([[1.0058, -0.0058], [-0.9829, 0.0056]])
Dh = np.array([[-1], [1]])

M = np.array([[0, 1], [1, 0]])
ni, n = M.shape[0] / 2, Ah.shape[0]
rho = 0.5

P = cvx.Semidef(n)
lamda = cvx.Variable()

Mblk = M*lamda*np.eye(2)
ABh = cvx.hstack(Ah, Bh)
CDh = cvx.hstack(Ch, Dh)
zeros = np.zeros((n,1))
constraints = [lamda[-1] == 1,
               ABh.T*P*ABh - rho**2*cvx.bmat([[P,zeros],[zeros.T, 0]]) +
               CDh.T*Mblk*CDh << 0]

prob = cvx.Problem(cvx.Minimize(1), constraints)
prob.solve()
feas = prob.status is cvx.OPTIMAL

Я удалил функцию kron, потому что она здесь ничего не делала, а CVXPY в настоящее время не поддерживает продукты Kronecker с переменной левой частью. Я могу добавить это, если вам это нужно.