有一个Python代码可以解决最优控制问题，但我没有得到解决方案。我不知道哪里出了问题

Question

目标是最小化最后的隔间

M_u

和

F_u

并最大化

M_i

和

F_i

。我的职能目标是：

在乳胶中：

$J(u)=M_u^2+F_u^2+ \int_0^T(Mi+Fi)u^2dt$

控制的颂歌系统如下。其中

u

是控件

在乳胶中：

\begin{align}
\alpha =& \left( 1-\frac{A}{K_a}\right)\\
\delta =& (\gamma + \mu_a )A_u \\
\sigma =& M_uF_u+(1-\rho_1) M_iF_i +(1-\rho_2)M_uF_i+(1-\rho_3\lambda_m) M_iF_u \\
\frac{dA_u}{dt}=&\phi\left( \frac{\sigma}{M_u+M_i}\right)\alpha-\delta \\
\frac{dA_i}{dt}=&\phi\left( \frac{\rho_1M_iF_i+\rho_2M_uF_i+ \rho_3\lambda_m M_iF_u}{M_u+ M_i}\right) \left( 1-\frac{A}{K_a}\right) - \delta \\
\frac{dF_u}{dt}=&b_f\gamma A_u -\beta_m \frac{M_i}{M_u+ M_i} F_u- \mu_f F_u \\
\frac{dF_i}{dt}=&b_f\gamma A_i +\beta_m \frac{M_i}{M_u+ M_i} F_u- \mu_f F_i+uF_i \\
\frac{dM_u}{dt}=&b_m\gamma A_u -\beta_f \frac{F_i}{F_u+F_i} M_u- \mu_m M_u \\
\frac{dM_i}{dt}=&b_m\gamma A_i +\beta_f \frac{F_i}{F_u+F_i} M_u- \mu_m M_i +uM_i.
\end{align}

这是Python代码：

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)

Au = m.Var(10000)
Ai = m.Var(0)
Fu = m.Var(15000)
Fi = m.Var(200)
Mu = m.Var(15000)
Mi = m.Var(200)
M = m.Var(1,lb=1)
F = 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)


Ka = 20000
phi =5
bf = 0.5
bm = 0.5
ro1= 0.76
ro2 = 0.45
ro3 = 0.37
lm =0.51
gamma =0.2
betam = 0.035
betaf = 0.032
ua = 0.025

um = 0.12
uf = 0.1

# system
m.Equation(M==Mu+Mi)
m.Equation(F==Fu+Fi)
m.Equation(Au.dt() == phi*(Mu*Fu+(1-ro1)*Mi*Fi+(1-ro2)*Mu*Fi+(1-ro3*lm)*Mi*Fu)*(1-(Au+Ai)/Ka)/M-(gamma+ua)*Au)
m.Equation(Ai.dt() == phi*(ro1*Mi*Fi+ro2*Mu*Fi+ro3*lm*Mi*Fu)*(1-(Au+Ai)/Ka)/M-(gamma+ua)*Ai)
m.Equation(Fu.dt() == bf*gamma*Au-betam*Mi*Fu/M-uf*Fu)
m.Equation(Fi.dt() == bf*gamma*Ai+betam*Mi*Fu/M-uf*Fi+u*Fi)
m.Equation(Mu.dt() == bm*gamma*Au-betaf*Fi*Mu/F-um*Mu)
m.Equation(Mi.dt() == bm*gamma*Ai+betaf*Fi*Mu/F-um*Mi+u*Mi)
#m.Equation(final*(y-0.0)==0)
m.Minimize(final*((Fu-0.0)**2+(Mu-0.0)**2))

# cost function
m.Minimize(m.integral((Mi+Fi)*u**2)*final)
           
m.options.IMODE = 6
m.options.MAX_ITER = 2000
m.options.NODES = 2
m.options.MV_TYPE = 1
m.options.SOLVER = 2

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.plot(m.time,Au.value+Fu.value+Mu.value,lw=2,label=r'$y$')
plt.plot(m.time,Ai.value+Fi.value+Mi.value,lw=2,label=r'$x$')
plt.grid()
plt.plot(m.time,u.value,'g',lw=2,label=r'$u$')
plt.grid()
plt.xlabel('Time')

有什么建议可以找到解决方案吗？