将半定程序从CVX转换为CVXPY [英] Convert a semidefinite program from CVX to CVXPY
问题描述
我想将以下 SDP(它只是验证约束的可行性)从 CVX (MATLAB) 转换为 CVXPY (Python):
I want to convert the following SDP — which just verifies the feasibility of the constraints — from CVX (MATLAB) to 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代码,
Below is my Python code,
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时,它会显示
There are several errors when I run the program. 1. When I print Mblk, it shows
回溯(最近通话最近一次):
Traceback (most recent call last):
文件"/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py",第2820行,在run_code中
File "/usr/lib/python2.7/dist-packages/IPython/core/interactiveshell.py", line 2820, in run_code
Out [1]:self.user_global_ns,self.user_ns中的exec code_obj
Out[1]: exec code_obj in self.user_global_ns, self.user_ns
文件",位于
兆位
文件"/usr/lib/python2.7/dist-packages/IPython/core/displayhook.py",第247行,在致电
File "/usr/lib/python2.7/dist-packages/IPython/core/displayhook.py", line 247, in call
format_dict,md_dict = self.compute_format_data(结果)
format_dict, md_dict = self.compute_format_data(result)
文件"/usr/lib/python2.7/dist-packages/IPython/core/displayhook.py",第157行,在compute_format_data中
File "/usr/lib/python2.7/dist-packages/IPython/core/displayhook.py", line 157, in compute_format_data
返回self.shell.display_formatter.format(结果)
return self.shell.display_formatter.format(result)
文件"/usr/lib/python2.7/dist-packages/IPython/core/formatters.py",第152行,格式为
File "/usr/lib/python2.7/dist-packages/IPython/core/formatters.py", line 152, in format
数据=格式化程序(obj)
data = formatter(obj)
文件/usr/lib/python2.7/dist-packages/IPython/core/formatters.py",第481行,在致电
File "/usr/lib/python2.7/dist-packages/IPython/core/formatters.py", line 481, in call
printer.pretty(obj)
printer.pretty(obj)
文件"/usr/lib/python2.7/dist-packages/IPython/lib/pretty.py",行362,漂亮
File "/usr/lib/python2.7/dist-packages/IPython/lib/pretty.py", line 362, in pretty
返回_default_pprint(obj,self,cycle)
return _default_pprint(obj, self, cycle)
文件"/usr/lib/python2.7/dist-packages/IPython/lib/pretty.py",行482,在_default_pprint
File "/usr/lib/python2.7/dist-packages/IPython/lib/pretty.py", line 482, in _default_pprint
p.text(repr(obj))
p.text(repr(obj))
文件"/usr/lib/python2.7/dist-packages/numpy/core/numeric.py",行1553,在array_repr中
File "/usr/lib/python2.7/dist-packages/numpy/core/numeric.py", line 1553, in array_repr
',',"array(")
', ', "array(")
文件"/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py",行454,在array2string中
File "/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py", line 454, in array2string
分隔符,前缀,formatter = formatter)
separator, prefix, formatter=formatter)
文件"/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py",行256,以_array2string
File "/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py", line 256, in _array2string
'int':IntegerFormat(数据),
'int' : IntegerFormat(data),
文件"/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py",行641,在 init
File "/usr/lib/python2.7/dist-packages/numpy/core/arrayprint.py", line 641, in init
max_str_len = max(len(str(max(maximum.reduce(data))),
max_str_len = max(len(str(maximum.reduce(data))),
文件"/usr/local/lib/python2.7/dist-packages/cvxpy/constraints/leq_constraint.py",第67行,非零
File "/usr/local/lib/python2.7/dist-packages/cvxpy/constraints/leq_constraint.py", line 67, in nonzero
引发异常(无法评估约束的真值.")
Raise Exception("Cannot evaluate the truth value of a constraint.")
异常:无法评估约束的真值.
Exception: Cannot evaluate the truth value of a constraint.
当我走到这一行时,
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 (most recent call last): File
".../sdp.py",第22行,位于
".../sdp.py", line 22, in
np.dot(ABh.T, np.dot(P, ABh)) - rho**2*np.linalg.block_diag(P, 0) +
ValueError:设置具有序列的数组元素.
ValueError: setting an array element with a sequence.
如何解决这些问题?
推荐答案
您的代码的最大问题是您不能在CVXPY对象上使用NumPy函数.您需要使用等效的CVXPY函数.这是您的代码的有效版本:
The big issue with your code is that you can't use NumPy functions on CVXPY objects. You need to use the equivalent CVXPY functions. Here's a working version of your code:
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产品.我可以根据需要添加它.
I removed the kron function because it wasn't doing anything here and CVXPY doesn't currently support Kronecker products with a variable left-hand side. I can add it if you need it.
这篇关于将半定程序从CVX转换为CVXPY的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!