GEKKO Multiphase Optimization: No solution, -1 degree of freedom

Question

I am learning GEKKO right now and tried to code up a simple 1D control example where a block with dynamics

dx/dt = v
dv/dt = u

with: -5 <= u <= 5

where x is the position, v is the velocity and u is the applied actuation has to move from position 0 to 1 in 1 time unit with initial and final velocity being 0. The objective is:

min abs(u)

It worked with a single phase and the result is this Bang-Off-Bang control.

I now wanted to try to implement this as a multiphase trajectory with 3 linked phases:

1. u = 5
2. u = 0
3. u = -5

where the time for each phase is the control variable.

The problem is that the solver does not find a solution. The degrees of freedom are -1 where they should be 1 in theory (the overall final time) with the given objective or 0 if the objective is replaced by an equality constraint (which results in the solver stating even less degrees of freedom, about -100).

Here is my code:

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO(remote=True)  # Model

n = 3  # Number of phases: Bang-Off-Bang
max_u = 5.0  # Control can either be max_u, 0 or -max_u

# Options
m.options.NODES = 4  # Two intermediate nodes between collocation points
m.options.SOLVER = 1  # APOPT
m.options.IMODE = 6  # MPC Direct Collocation
m.options.MAX_ITER = 500  # To prevent endless loops
m.options.MV_TYPE = 0  # MVs are constant between endpoints
m.options.DIAGLEVEL = 1  # Show some diagnostics

# Time from 0 to 1 in 31 collocation points for every phase each
m.time = np.linspace(0, 1, 31)

tf = [m.FV(value=0.3) for i in range(n)]
for i in range(n):
    tf[i].STATUS = 1  # make final time controllable

# Collocation variables x and v
x = [m.Var(value=0, fixed_initial=False) for i in range(n)]
v = [m.Var(value=0, fixed_initial=False) for i in range(n)]
u = [max_u, 0, -max_u]  # For the three phases

# Fix start- and endpoint (Inbetween points can be fixed with m.fix())
m.fix_initial(x[0], val=0)
m.fix_initial(v[0], val=0)
m.fix_final(x[n-1], val=1)
m.fix_final(v[n-1], val=0)

# Differential equations describing the system scaled by tf
for i in range(n):
    m.Equation(x[i].dt() == v[i]*tf[i])
    m.Equation(v[i].dt() == u[i]*tf[i])

# Connect phases at endpoints
for i in range(n-1):
    m.Connection(x[i+1], x[i], 1, 'end', 1, 'end')
    m.Connection(x[i+1],'calculated', pos1=1, node1=1)
    m.Connection(v[i+1], v[i], 1, 'end', 1, 'end')
    m.Connection(v[i+1],'calculated', pos1=1, node1=1)

# Make final time equal to 1
m.Minimize((m.sum(tf)-1)**2)
#m.Equation(m.sum(tf) == 1)  # Could be used instead of objective

m.open_folder()

# Run optimization with Diagnostics
m.solve(disp=True)

I hope someone can help me, thank you in advance!

Answer 1

Soft terminal constraints typically do better for converging to a solution. Here is a replacement of the hard terminal constraint:

f = np.zeros(31); f[-1]=1; final=m.Param(f)
m.Minimize(final*(x[n-1]-1)**2)
m.Minimize(final*v[n-1]**2)

#m.fix_final(x[n-1], val=1)
#m.fix_final(v[n-1], val=0)

This gives an optimal solution with the three linked phases.

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

m = GEKKO(remote=True)  # Model

n = 3  # Number of phases: Bang-Off-Bang
max_u = 5.0  # Control can either be max_u, 0 or -max_u

# Options
m.options.NODES = 4  # Two intermediate nodes between collocation points
m.options.SOLVER = 1  # APOPT
m.options.IMODE = 6  # MPC Direct Collocation
m.options.MAX_ITER = 500  # To prevent endless loops
m.options.MV_TYPE = 0  # MVs are constant between endpoints
m.options.DIAGLEVEL = 0  # Show some diagnostics

# Time from 0 to 1 in 31 collocation points for every phase each
m.time = np.linspace(0, 1, 31)

tf = [m.FV(value=0.3) for i in range(n)]
for i in range(n):
    tf[i].STATUS = 1  # make final time controllable

# Collocation variables x and v
x = [m.Var(value=0, fixed_initial=False) for i in range(n)]
v = [m.Var(value=0, fixed_initial=False) for i in range(n)]
u = [max_u, 0, -max_u]  # For the three phases

# Fix start- and endpoint (Inbetween points can be fixed with m.fix())
m.fix_initial(x[0], val=0)
m.fix_initial(v[0], val=0)

f = np.zeros(31); f[-1]=1; final=m.Param(f)
m.Minimize(final*(x[n-1]-1)**2)
m.Minimize(final*v[n-1]**2)

#m.fix_final(x[n-1], val=1)
#m.fix_final(v[n-1], val=0)

# Differential equations describing the system scaled by tf
for i in range(n):
    m.Equation(x[i].dt() == v[i]*tf[i])
    m.Equation(v[i].dt() == u[i]*tf[i])

# Connect phases at endpoints
for i in range(n-1):
    m.Connection(x[i+1], x[i], 1, 'end', 1, 'end')
    m.Connection(x[i+1],'calculated', pos1=1, node1=1)
    m.Connection(v[i+1], v[i], 1, 'end', 1, 'end')
    m.Connection(v[i+1],'calculated', pos1=1, node1=1)

# Make final time equal to 1
m.Minimize((m.sum(tf)-1)**2)
#m.Equation(m.sum(tf) == 1)  # Could be used instead of objective

#m.open_folder()

# Run optimization with Diagnostics
m.solve(disp=True)

# Generate plot
t = [m.time*tf[i].value[0] for i in range(3)]
t[1] += t[0][-1]
t[2] += t[1][-1]
for i in range(3):
    plt.plot(t[i],x[i].value)
    plt.plot(t[i],v[i].value)
plt.xlabel('Time')
plt.show()    

The missing 2 degrees of freedom (-1 DOF) is because the derivatives at the endpoint are also fixed when m.fix_final(x[n-1], val=1) and m.fix_final(v[n-1], val=0) are used. This is a known issue with Gekko and how the variables are defined. There is an alternative way (create a fixed FV that connects to the endpoint) if you do need a hard constraint.