3

I'm Arthur. I'm very new to gekko, thus I might be missing something very obvious here.

I have been trying to find the solution to an optimal control problem, namely trajectory optimization of a spacecraft rendezvous problem (Chaser -> Target), under certain final constraints to the relative final distance and velocity and certain distance along their trip. In order to do this, I tried using soft constraints to the final path. As a objective function, I use a norm 1 of the force generated from Thrusters.

The computation keeps going until the maximum no. of iterations (2000) is reached and cancels. Is there maybe a way to discretize the search space of this problem to make it solvable in acceptable time trading off computation time for accuracy?

import matplotlib.pyplot as plt
import numpy as np

from gekko import GEKKO

# create GEKKO model
m = GEKKO()
print(m.path)
nt = 501
m.time = np.linspace(0,500,nt)

# Variables

# Initial Position and Velocity - Chaser
r1_c = -1182.339411   # [km] 
r2_c = 6816.939420    # [km] 
r3_c = 904.891745     # [km] 
v1_c = 0.1175776      # [km/s] 
v2_c = -0.963776      # [km/s] 
v3_c = 7.494102       # [km/s] 

# Initial Position and Velocity - Target
r1_t = -1182.959348    # [km] 
r2_t = 6817.396210     # [km] 
r3_t = 904.495486      # [km] 
v1_t = 0.1175776       # [km/s] 
v2_t = -0.963776       # [km/s] 
v3_t = 7.494102        # [km/s] 

# initial values
x =  m.Var(value = r1_c)
y =  m.Var(value = r2_c)
z =  m.Var(value = r3_c)
vx =  m.Var(value = v1_c)
vy =  m.Var(value = v2_c)
vz =  m.Var(value = v3_c)


# initial values
x2 =  m.Var(value = r1_t)
y2 =  m.Var(value = r2_t)
z2 =  m.Var(value = r3_t)
vx2 =  m.Var(value = v1_t)
vy2 =  m.Var(value = v2_t)
vz2 =  m.Var(value = v3_t)

# Cost Initial - Integrated thrust (fuel usage)
J =  m.Var(value = 0)

# Manipulated Variable
U = m.MV(value=0, lb=-0.001, ub=0.001)

U.STATUS = 1

# Mask time
p = np.zeros(nt) # mask final time point
p[-1] = 500
final = m.Param(value=p)

# Parameters
RT =  m.Const(6378.139);            # [km] 
mu_t = m.Const(3.9860064*10**(5))   # [km^3/s^2] 


# Equations

# Define intermediate quantities - Chaser
Rc = m.Intermediate((x**2 + y**2 + z**2)**0.5)
v = m.Intermediate((vx**2 + vy**2 + vz**2)**0.5)

ax = m.Intermediate(x * -mu_t / Rc**3 + U *vx / v )
ay = m.Intermediate(y * -mu_t / Rc**3 + U *vy / v )
az = m.Intermediate(z * -mu_t / Rc**3 + U *vz / v )


# Define intermediate quantities - Target
Rt = m.Intermediate((x2**2 + y2**2 + z2**2)**0.5)
v2 = m.Intermediate((vx2**2 + vy2**2 + vz2**2)**0.5)

ax2 = m.Intermediate(x2 * -mu_t / Rt**3 )
ay2 = m.Intermediate(y2 * -mu_t / Rt**3 )
az2 = m.Intermediate(z2 * -mu_t / Rt**3 )


# Governing equations
m.Equations((vx.dt() == ax, vy.dt() == ay, vz.dt() == az))
m.Equations((x.dt() == vx, y.dt() == vy, z.dt() == vz))

m.Equations((vx2.dt() == ax2, vy2.dt() == ay2, vz2.dt() == az2))
m.Equations((x2.dt() ==  vx2, y2.dt() ==  vy2, z2.dt() == vz2))

# Equation relating thrust to fuel usage
m.Equation(J.dt() == m.abs2(U))


# Path Constraints

# specify endpoint conditions

m.fix(x, pos=len(m.time)-1, val=x2)
m.fix(y, pos=len(m.time)-1, val=y2)
m.fix(z, pos=len(m.time)-1, val=z2)
m.fix(vx, pos=len(m.time)-1, val=vx2)
m.fix(vy, pos=len(m.time)-1, val=vy2)
m.fix(vz, pos=len(m.time)-1, val=vz2)

# Constraints 

m.Equations((m.abs2(Rc - RT) > 200, m.abs2(Rc - RT) < 1000, m.abs2(Rc*final - Rt*final) == 0, m.abs2(v*final - v2*final) == 0))


# Objective, minimizing fuel usage
m.Obj(J * final)

# Set solver mode - Optimal Control
m.options.IMODE = 6

# Increase maximum number of allowed iterations
m.options.MAX_ITER = 2000

# Set number of nodes per time segment
m.options.NODES = 3

# Run solver and display progress
m.solve(disp=True)

---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :    3663.92820000000      sec
 Objective      :    185346862893.104     
 Unsuccessful with error code            0
 ---------------------------------------------------
 
 Creating file: infeasibilities.txt
 Use command apm_get(server,app,'infeasibilities.txt') to retrieve file
 @error: Solution Not Found
---------------------------------------------------------------------------

It does not find the solution and I do not find the file infeasibilities.txt.

Can you help me please?

I would like to find the reason from results (Solution Not Found). However someone can help me with structure of the problem becoming it feasible.

John Hedengren
  • 12,068
  • 1
  • 21
  • 25
Arthur Moreno
  • 31
  • 1
  • Here is some help on retrieving the infeasibilities.txt report: https://stackoverflow.com/questions/56671888/how-to-retrieve-the-infeasibilities-txt-from-the-gekko – TexasEngineer Jan 30 '23 at 16:20

1 Answers1

0

One infeasible set of constraints is that Rc-RT is set to >200 and <1000 everywhere, but it must also equal zero at the final time. The solver can't satisfy these constraints.

m.abs2(Rc - RT) > 200
m.abs2(Rc - RT) < 1000
m.abs2(Rc*final - Rt*final) == 0

I recommend replacing these constraints with something that is numerically simple with the definition of a new variable with upper and lower bounds and the **2 to avoid sqrt() and abs().

RcRT = m.Var(lb=200,ub=1000)
m.Equation(RcRT**2 == (Rc-RT)**2)

The MPCC m.abs2() does not work very well for most problems. m.abs3() works better, but can take more time to solve as a mixed integer form. Try to avoid the use of abs() in the model, when there is an alternative. Equations such as m.abs2(v*final - v2*final) == 0 will have no solution if the Chaser is not able to match the velocity of the Target. Try a different soft constraint (objective-based) that allows the problem to remain feasible and minimize the difference.

m.Minimize(final*(v-v2)**2)

Likewise, these hard constraints may be infeasible:

# hard constraints
m.fix(x, pos=len(m.time)-1, val=x2)
m.fix(y, pos=len(m.time)-1, val=y2)
m.fix(z, pos=len(m.time)-1, val=z2)
m.fix(vx, pos=len(m.time)-1, val=vx2)
m.fix(vy, pos=len(m.time)-1, val=vy2)
m.fix(vz, pos=len(m.time)-1, val=vz2)

Try soft constraints:

# soft constraints
m.Minimize(final*(v-v2)**2)
m.Minimize(final*(x-x2)**2)
m.Minimize(final*(y-y2)**2)
m.Minimize(final*(z-z2)**2)
m.Minimize(final*(vx-vx2)**2)
m.Minimize(final*(vy-vy2)**2)
m.Minimize(final*(vz-vz2)**2)

Here is the complete problem:

import matplotlib.pyplot as plt
import numpy as np

from gekko import GEKKO

# create GEKKO model
m = GEKKO()
print(m.path)
nt = 101
m.time = np.linspace(0,500,nt)

# Variables

# Initial Position and Velocity - Chaser
r1_c = -1182.339411   # [km] 
r2_c = 6816.939420    # [km] 
r3_c = 904.891745     # [km] 
v1_c = 0.1175776      # [km/s] 
v2_c = -0.963776      # [km/s] 
v3_c = 7.494102       # [km/s] 

# Initial Position and Velocity - Target
r1_t = -1182.959348    # [km] 
r2_t = 6817.396210     # [km] 
r3_t = 904.495486      # [km] 
v1_t = 0.1175776       # [km/s] 
v2_t = -0.963776       # [km/s] 
v3_t = 7.494102        # [km/s] 

# initial values
x =  m.Var(value = r1_c)
y =  m.Var(value = r2_c)
z =  m.Var(value = r3_c)
vx =  m.Var(value = v1_c)
vy =  m.Var(value = v2_c)
vz =  m.Var(value = v3_c)


# initial values
x2 =  m.Var(value = r1_t)
y2 =  m.Var(value = r2_t)
z2 =  m.Var(value = r3_t)
vx2 =  m.Var(value = v1_t)
vy2 =  m.Var(value = v2_t)
vz2 =  m.Var(value = v3_t)

# Cost Initial - Integrated thrust (fuel usage)
J =  m.Var(value = 0)

# Manipulated Variable
U = m.MV(value=0, lb=-0.001, ub=0.001)
U.STATUS = 1

# Mask time
p = np.zeros(nt) # mask final time point
p[-1] = 500
final = m.Param(value=p)

# Parameters
RT =  m.Const(6378.139);            # [km] 
mu_t = m.Const(3.9860064*10**(5))   # [km^3/s^2] 

# Equations

# Define intermediate quantities - Chaser
Rc = m.Intermediate((x**2 + y**2 + z**2)**0.5)
v = m.Intermediate((vx**2 + vy**2 + vz**2)**0.5)

ax = m.Intermediate(x * -mu_t / Rc**3 + U *vx / v )
ay = m.Intermediate(y * -mu_t / Rc**3 + U *vy / v )
az = m.Intermediate(z * -mu_t / Rc**3 + U *vz / v )

# Define intermediate quantities - Target
Rt = m.Intermediate((x2**2 + y2**2 + z2**2)**0.5)
v2 = m.Intermediate((vx2**2 + vy2**2 + vz2**2)**0.5)

ax2 = m.Intermediate(x2 * -mu_t / Rt**3 )
ay2 = m.Intermediate(y2 * -mu_t / Rt**3 )
az2 = m.Intermediate(z2 * -mu_t / Rt**3 )


# Governing equations
m.Equations((vx.dt() == ax, vy.dt() == ay, vz.dt() == az))
m.Equations((x.dt() == vx, y.dt() == vy, z.dt() == vz))

m.Equations((vx2.dt() == ax2, vy2.dt() == ay2, vz2.dt() == az2))
m.Equations((x2.dt() ==  vx2, y2.dt() ==  vy2, z2.dt() == vz2))

# Equation relating thrust to fuel usage
#m.Equation(J.dt() == m.abs2(U))
m.Equation(J.dt()**2 == U**2)


# Path Constraints

# specify endpoint conditions
# soft constraints
m.Minimize(final*(v-v2)**2)
m.Minimize(final*(x-x2)**2)
m.Minimize(final*(y-y2)**2)
m.Minimize(final*(z-z2)**2)
m.Minimize(final*(vx-vx2)**2)
m.Minimize(final*(vy-vy2)**2)
m.Minimize(final*(vz-vz2)**2)    

# Constraints 
# m.abs2(Rc - RT) > 200
# m.abs2(Rc - RT) < 1000
# m.abs2(Rc*final - Rt*final) == 0
RcRT = m.Var(lb=200,ub=1000)
m.Equation(RcRT**2 == (Rc-RT)**2)

# Objective, minimizing fuel usage
m.Minimize(J * final)

# Set solver mode - Optimal Control
m.options.IMODE = 6

# Increase maximum number of allowed iterations
m.options.MAX_ITER = 2000

# Set number of nodes per time segment
m.options.NODES = 3

# Run solver and display progress
m.solve(disp=True)

IPOPT finds an optimal solution after 1375 iterations with a reduced the number of time points to 101 for testing. You can try to increase the number of time points and resolve.

1375 -1.3996401e+02 2.04e-09 8.64e-07 -10.5 1.61e-06    -  1.00e+00 1.00e+00h  1
Number of Iterations....: 1375
                                   (scaled)                 (unscaled)
Objective...............:  -4.4944986703206206e-03   -1.3996401074391898e+02
Dual infeasibility......:   8.6432001202008113e-07    2.6915948656834957e-02
Constraint violation....:   7.7449958484773991e-10    2.0372681319713593e-09
Complementarity.........:   3.3628245954861400e-11    1.0472233998435268e-06
Overall NLP error.......:   8.6432001202008113e-07    2.6915948656834957e-02

Number of objective function evaluations             = 1974
Number of objective gradient evaluations             = 1149
Number of equality constraint evaluations            = 1982
Number of inequality constraint evaluations          = 1982
Number of equality constraint Jacobian evaluations   = 1389
Number of inequality constraint Jacobian evaluations = 1389
Number of Lagrangian Hessian evaluations             = 1375
Total CPU secs in IPOPT (w/o function evaluations)   =     83.311
Total CPU secs in NLP function evaluations           =    112.375

EXIT: Optimal Solution Found.
 The solution was found.
 The final value of the objective function is   -139.964010743919     
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :    199.320599999999      sec
 Objective      :   -139.964010743919     
 Successful solution
 ---------------------------------------------------
John Hedengren
  • 12,068
  • 1
  • 21
  • 25