0

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

よろしくお願いいたします。

4

1 に答える 1