Hello everyone. I would like to solve an optimal control problem by minimizing an integral but I get solution not found. here is the code [closed]

Question

I would like to solve an optimal control problem by minimizing an integral, but I get solution not found. Here is the code.

J(u)=\int{0}^{500}(69732*u*x+u**2)dt with the constraint y(500)=0

import numpy as np
from gekko import GEKKO
m = GEKKO(remote=False)
nt = 101; m.time = np.linspace(0,500,nt)

y = m.Var(15000)
x = m.Var(150)
u = m.MV(value=0.06,lb=0.02, ub=0.06)
u.STATUS=1
p = np.zeros(nt); p[-1] = 1.0; final = m.Param(value=p)

# system
m.Equation(y.dt()==3.5*(y*y+0.55*x*y+(1-0.37*0.53)*x*y)*(1-(x+y)/20000)/(x+y)-0.035*x*y/(x+y)-0.1*y)
m.Equation(x.dt()==3.5*(x*x+0.45*x*y+0.37*0.53*x*y)*(1-(x+y)/20000)/(x+y)+0.035*x*y/(x+y)-0.1*x+u*x)
m.Equation(final*(y-0.0)==0)

# cost function
m.Minimize(m.integral(69732*u*x+u**2)*final)
           
m.options.IMODE = 6
m.options.NODES = 2
m.options.MV_TYPE = 1
m.options.SOLVER = 3
m.solve()
print('Objective: ' + str(m.options.OBJFCNVAL))

import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(m.time,y.value,'r:',LineWidth=2,label=r'$y$')
plt.plot(m.time,x.value,'k:',LineWidth=2,label=r'$x$')
plt.xlabel('Time')
plt.grid(which="both")
plt.legend(loc='best')
plt.figure(2)
plt.plot(m.time,u.value,'b--',LineWidth=2,label=r'$u$')
plt.xlabel('Time')
plt.grid(which="both")
plt.show()```

Answer 1

The divide-by-zero with x+y is a potential problem as well as the inability to reach the terminal constraint y(500)=0. Here is a modified problem that solves successfully with a new variable z and a softened terminal constraint minimize(y(500)**2).

import numpy as np
from gekko import GEKKO
m = GEKKO(remote=True)

m.time = np.concatenate((np.array([0.0,0.01,0.02,0.05,0.1,0.2,0.5]),\
                        np.linspace(1,500,500)))
nt = len(m.time)

y = m.Var(15000)
x = m.Var(150)
z = m.Var(1,lb=1)
u = m.MV(value=0.06,lb=0.02, ub=0.06)
u.STATUS=1
p = np.zeros(nt); p[-1] = 1.0; final = m.Param(value=p)

# system
m.Equation(z==x+y)
m.Equation(y.dt()==3.5*(y*y+0.55*x*y+(1-0.37*0.53)*x*y)*(1-z/20000)/z-0.035*x*y/z-0.1*y)
m.Equation(x.dt()==3.5*(x*x+0.45*x*y+0.37*0.53*x*y)*(1-z/20000)/z+0.035*x*y/z-0.1*x+u*x)
#m.Equation(final*(y-0.0)==0)
m.Minimize(final*(y-0.0)**2)

# cost function
m.Minimize(m.integral(69732*u*x+u**2)*final)
           
m.options.IMODE = 6
m.options.MAX_ITER = 1000
m.options.NODES = 2
m.options.MV_TYPE = 1
m.options.SOLVER = 3

m.options.COLDSTART = 1
m.solve()

m.options.COLDSTART = 0
m.options.TIME_SHIFT = 0
m.solve()

print('Objective: ' + str(m.options.OBJFCNVAL))

import matplotlib.pyplot as plt
plt.figure(1)
plt.subplot(3,1,1)
plt.plot(m.time,y.value,'r:',lw=2,label=r'$y$')
plt.grid(); plt.legend()
plt.subplot(3,1,2)
plt.plot(m.time,x.value,'k:',lw=2,label=r'$x$')
plt.grid(); plt.legend()
plt.subplot(3,1,3)
plt.plot(m.time,u.value,'b--',lw=2,label=r'$u$')
plt.xlabel('Time')
plt.grid(); plt.legend()
plt.savefig('results.png',dpi=300)
plt.show()