How to deal with the transversality condition (adjoint variables Lamda) in Python Gekko Optimization?

Here are a few suggestions:

1. La1 - La5 need to be variables instead of fixed parameters.

La1 = m.Var()
La2 = m.Var()
La3 = m.Var()
La4 = m.Var()
La5 = m.Var()

You could also define them on one line with m.Array().

La1,La2,La3,La4,La5 = m.Array(m.Var,5)

2. The constants a0 - e4 need to be m.Param(), m.Const(), or just python floating point numbers. If they are parameters that need to be calculated, set them as a0=m.FV(lb=0,ub=1) types with lower and upper bounds and a0.STATUS=1 to make it a degree of freedom for the solver.

a0 = 1e-4
b0 = 0.1
c0 = 0.005
f0 = 1
a1 = 1e-4
b1 = 1e-2
c1 = 0.005
f1 = 1
b2 = 0.02
d2 = 0.1
f2 = 1
d3 = 0.1
b4 = 1e-2
c4 = 1e-2
e4 = 1e-7

3. Minimize only the last value of the integral. Fix this by adding parenthesis.

m.Minimize(tf*(M+0.5*m.integral(u**2)))

4. Switch to IMODE=4 or IMODE=7 (faster) for simulation. Turn off the display when running IMODE=7 to not waste time with the large printout.

m.options.IMODE=7
m.solve(disp=False)

5. Equations such as m.Equation(La3*tf==1) are infeasible if tf=0 because it gives an equation 0=1. Instead, use m.Equation(tf*(La3-1)==0) so that the equation is only active at the end. Before those equations can be added, additional degrees of freedom need to be added. Another way to write the equations is to use an objective function form such as m.Minimize(tf*(La3-1)**2). This is often easier for the solver to find a solution and won't return an error if the solver can't reach the final value.

from gekko import GEKKO
import numpy as np

m = GEKKO()

nt = 501
m.time = np.linspace(0,500,nt)

# Variables
u = m.MV(value=0,lb=0,ub=1)
u.STATUS = 1
H = m.Var(value=0.1,lb=0)
C = m.Var(value=0.1,lb=0)
M = m.Var(value=0.8,lb=0)
D = m.Var(value=0.1,lb=0)
I = m.Var(value=0.1,lb=0)

a0 = 1e-4
b0 = 0.1
c0 = 0.005
f0 = 1
a1 = 1e-4
b1 = 1e-2
c1 = 0.005
f1 = 1
b2 = 0.02
d2 = 0.1
f2 = 1
d3 = 0.1
b4 = 1e-2
c4 = 1e-2
e4 = 1e-7

p = np.zeros(nt)
p[-1] = 1
tf = m.Param(value=p)
La1,La2,La3,La4,La5 = m.Array(m.Var,5)

# Control Variables Equations
m.Equation(H.dt()==a0-c0*H+b0*D*H*(1-H/f0))
m.Equation(C.dt()==a1-c1*C+b1*I*(M+D)*C*(1-C/f1))
m.Equation(M.dt()==b2*M*(1-M/f2)-d2*M*C)
m.Equation(D.dt()==-d3*D*C+u)
m.Equation(I.dt()==b4*D*H-e4*I*C+c4*I)

# Adjoint Variables Equations
m.Equation(La1.dt()==(-c0+b0*D-2*D*b0*H/f0)*La1+b4*D*La5)
m.Equation(La2.dt()==(-c1+b1*I*(M+D)*(1-2*C/f1))*La2-d2*M*La3-d3*D*La4-e4*I*La5)
m.Equation(La3.dt()==-1+b1*I*C*(1-C/f1)*La2+(b2*(1-2*M/f1)-d2*C)*La3)
m.Equation(La4.dt()==b1*H*(1-H/f0)*La1+b1*I*C*(1-C/f1)*La2-d3*C*La4+b4*H*La5)
m.Equation(La5.dt()==b1*(M+D)*C*(1-C/f1)*La2+(-e4*C-c4)*La5)
#m.Equation(La1*tf==0)
#m.Equation(La2*tf==0)
#m.Equation(La3*tf==1)
#m.Equation(La4*tf==0)
#m.Equation(La5*tf==0)           

# Objective Function
m.Minimize(tf*(M+0.5*m.integral(u**2)))

m.options.IMODE=7
m.solve(disp=False)

import matplotlib.pyplot as plt
plt.figure(figsize=(7,4))
plt.plot(m.time,H,label='H')
plt.plot(m.time,C,label='C')
plt.plot(m.time,M,label='M')
plt.plot(m.time,D,label='D')
plt.plot(m.time,I,label='I')
plt.legend()
plt.savefig('result.png',dpi=300)
plt.show()

Don't forget to include the [gekko] tag to get noticed earlier when asking a new question.