我想使用 CVXPY 在离散时间解决最小燃料最优控制问题。在连续时间线性系统上使用零阶保持，可以将问题转化为具有凸控制约束的线性程序。我已经使用 Yalmip 和 CVX 等 Matlab 环境完成了这个问题的基本公式，但我无法在 CVXPY 中使用它。我遇到的问题是，即使问题似乎可以编译并完全解决，但输出值显然不满足边界条件，并且在绘制输出时，结果为空。
我附上了代码。问题应该是双积分器的最小燃料控制；我想要一个严格正向和负向的控制信号，每个信号的取值都在 0 和 umax 之间。
这是定义双积分器矩阵的类：
import numpy as np
import control as ctrl
class DoubleIntegrator:
def __init__(self, nu):
self.A = np.matrix([[0,1],[0,0]])
# Note, this matrix doesn't really do anything right now
self.C = np.matrix([[1,0],[0,1]])
try:
if nu == 1:
self.B = np.matrix([[0],[1]])
self.D = np.matrix([[0],[1]])
return
elif nu == 2:
self.B = np.matrix([[0,0],[1,-1]])
self.D = np.matrix([[0,0],[1,-1]])
return
elif nu != 1 or nu != 2:
raise InputError("ERROR: Input 'nu' must be either nu = 1 or nu = 2")
except InputError as me:
print me.message
return
def makeStateSpaceSystem(self):
self.sysc = ctrl.matlab.ss(self.A, self.B, self.C, self.D)
def makeDiscreteSystem(self, dt, method):
self.sysd = ctrl.matlab.c2d(self.sysc, dt, method)
class Error(Exception):
pass
class InputError(Error):
def __init__(self, message):
self.message = message
if __name__ == '__main__':
db = DoubleIntegrator(2)
db.makeStateSpaceSystem()
db.makeDiscreteSystem(0.1, 'zoh')
print db.sysd.A[0,1]
这是主要代码：
from DoubleIntegrator import *
import numpy as np
import cvxpy as cvx
import control as ctrl
import matplotlib.pyplot as plt
db = DoubleIntegrator(2)
db.makeStateSpaceSystem()
db.makeDiscreteSystem(0.1,'zoh')
A = db.sysd.A
B = db.sysd.B
Nsim = 20
x = cvx.Variable(2, Nsim+1)
u = cvx.Variable(2, Nsim)
x0 = [1.0,0.0]
xf = [0.0,0.0]
umax = 1.0
states = []
cost = 0
for t in range(Nsim):
cost = cost + u[0,t] + u[1,t]
constr = [x[:,t+1] == A*x[:,t] + B*u[:,t],
0 <= u[0,t], u[0,t] <= umax,
0 <= u[1,t], u[1,t] <= umax]
states.append( cvx.Problem( cvx.Minimize(cost), constr ) )
prob = sum(states)
prob.constraints += [x[0,Nsim] == xf[0], x[1,Nsim] == xf[1], x[0,0] == x0[0], x[1,0] == x0[1]]
prob.solve(verbose = True)
print u.value
print x.value
f = plt.figure()
plt.subplot(411)
plt.plot(x[0,:].value)
plt.subplot(412)
plt.plot(x[1,:].value)
plt.subplot(413)
plt.plot(u[0,:].value)
plt.subplot(414)
plt.plot(u[1,:].value)
plt.show()
非常感谢您的帮助！