Minimum fuel control with CVXPY
<p>I want to solve a minimum fuel optimal control problem in discrete time using CVXPY. Using a zero-order hold on a continuous time, linear system, the problem can be turned into a linear program with convex control constraints. I have done the basic formulation of this problem using Matlab environments like Yalmip and CVX, but I cannot get this to work in CVXPY. The issue I am having is that even though the problem seems to compile and solve completely, the output values clearly do not satisfy the boundary conditions, and when the output is plotted, the results are empty.</p>

<p>I have attached the code. The problem is supposed to be minimum fuel control of the double integrator; I want a control signal in the strictly positive and negative directions, each taking values between 0 and umax.</p>

<p>Here's the class defining the double integrator matrices:</p>

<pre><code>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]
</code></pre>

<p>Here's the main code:</p>

<pre><code>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 &lt;= u[0,t], u[0,t] &lt;= umax,
                   0 &lt;= u[1,t], u[1,t] &lt;= 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()
</code></pre>

<p>Thank you very much in advance for your help!</p>