Как исправить ошибку «Решение не найдено» в коде оптимального управления Python GEKKO - PullRequest
Новая
       23

Как исправить ошибку «Решение не найдено» в коде оптимального управления Python GEKKO

0 голосов
/ 09 июля 2019

Я пытаюсь воспроизвести результат на рисунке 1 статьи «Иммунотерапия: подход теории оптимального управления» К. Рене Фистер и Дженнифер Хьюз Доннелли, 2005 г. Для этого я написал числовое решение для оптимального управленияиспользуя пакет Python для GEKKO.Я использовал те же начальные условия, контрольные границы, значения параметров и модельные уравнения, что и в статье.Однако, когда я запускаю код, я получаю следующую ошибку: 

Traceback (most recent call last):
  File "xxxxxx", line 45, in <module>
    m.solve(disp=False) #solve
  File "xxxx", line 783, in solve
    raise Exception(response)
Exception:  @error: Solution Not Found

Я ожидаю, что выходные данные программы предоставят две цифры: одну из динамики ODE и график решения оптимального управления.

Я пытался изменить код разными способами: изменив целевой функционал, количество временных шагов и изменив оптимальный режим управления, однако я каждый раз получаю одну и ту же ошибку.Ниже приведен код, который я использую: 

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

m = GEKKO()
nt = 1010
m.time = np.linspace(0,350,nt)

# Variables
X = m.Var(value=1)
Y = m.Var(value=1)
Z = m.Var(value=1)
OF = m.Var(value=0)
v = m.Var(value=0,lb=0,ub=1) #Control is initially 0 with a lower bound of 0 and an upper bound of 1

p = np.zeros(nt) #mark final time point
p[-1] = 1.0 #all zeros except the end, which is 1
final = m.Param(value=p) #final depends on integration limits

#Parameter Values
c = .000085
mu2 = .03
p1 = .1245
a = 1
r2 = .18
mu3 = 10
p2 = 5
g1 = 20000000 #2e7
g2 = 100000 #1e5
g3 = 1000 #1e3
b = 1*10**(-9)
s2 = 100000000
B = 100000000

# Equations
m.Equation(X.dt() == c*Y-mu2*X+(p1*X*Z)/(g1+Z))
m.Equation(Y.dt() == r2*Y*(1-b*Y)-(a*X*Y)/(g2+Y))
m.Equation(Z.dt() == (p2*X*Y)/(g3+Y)-mu3*Z+v*s2)
m.Equation(OF.dt() == X-Y+Z-B*v)

m.Obj(-OF*final)

m.options.IMODE = 6 #optimal control mode
m.solve(disp=False) #solve

plt.figure(figsize=(4,3)) #plot results
plt.subplot(2,1,1)
plt.plot(m.time,X.value,'k-',label=r'$S$')
plt.plot(m.time,Y.value,'b-',label=r'$R$')
plt.plot(m.time,Z.value,'g-',label=r'$E$')
plt.plot(m.time,OF.value,'r-',label=r'$OF$')
plt.legend()
plt.ylabel('CV')
plt.subplot(2,1,2)
plt.plot(m.time,v.value,'g--',label=r'$v$')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

Этот код был получен путем изменения примера кода GEKKO, который был предоставлен в этом видео Youtube .Любая помощь в решении этой проблемы будет принята с благодарностью!

...