Welcome toVigges Developer Community-Open, Learning,Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
802 views
in Technique[技术] by (71.8m points)

python - Gekko: Problem with the obtained solution

I am trying to solve the following Optimal Control Problem in GEKKO:

problem

I know a priori that the path for c is:

enter image description here

where the parameter values are: r = 0.33, i = 0.5, K(0) = 10 and T = 10.

I wrote the following code in Python:

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

m = GEKKO(remote=True)
nt = 101; m.time = np.linspace(0,10,nt)
r = 0.33
i = 0.5

# Variables
c = m.Var()
k = m.Var(value=10)
objective = m.Var()

rate = [-r*t/10 for t in range(0, 101)]
factor = m.exp(rate)

p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)
disc = m.Param(value=factor)


# Equations
m.Equation(k.dt() == i*k - c)
m.Equation(objective.dt() == disc*m.log(c))
# Objective Function
m.Maximize(final*objective)

m.options.IMODE = 6
m.solve()

plt.figure(1)
plt.plot(m.time,c.value,'k:',LineWidth=2,label=r'$C$')
plt.plot(m.time,k.value,'b-',LineWidth=2,label=r'$K$')

plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

The solved path for c and k is as shown below: enter image description here

This clearly is not right because c should be increasing with time from just eye-balling the solution given before hand.

I'm not sure where I am wrong.


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Answer

0 votes
by (71.8m points)

The optimal control problem as it is currently written is unbounded. The value of c will go to infinity to maximize the function. I set an upper bound of 100 on c and the solver went to that bound. I reformulated the model to reflect the current problem statement. Here are a few suggestions:

  1. Use the m.integral() function to make the model more readable.
  2. Initialize c at a value other than 0 (default). You may also want to set a lower bound with c>0.01 so that m.log(c) is not undefined if the solver tries a value <0.
  3. Only use Gekko functions inside Gekko equations such as with factor = m.exp(rate). Use factor = np.exp(rate) instead unless it is in a Gekko equation where it can be evaluated.
  4. Include a plot of the exact solution so that you can compare the exact and numerical solution.
  5. Use m.options.NODES=3 with c=m.MV() and c.STATUS=1 to increase the solution accuracy. The default is m.options.NODES=2 that isn't as accurate.
  6. You can free the initial condition with m.free_initial(c) to calculate the initial value of c.

Unbounded Solution

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

m = GEKKO(remote=True)
nt = 101; m.time = np.linspace(0,10,nt)
r = 0.33
i = 0.5

# Variables
c = m.MV(4,lb=0.01,ub=100); c.STATUS=1
#m.free_initial(c)
k = m.Var(value=10)
objective = m.Var(0)
t = m.Param(m.time)
m.Equation(objective==m.exp(-r*t)*m.log(c))

# just to include on the plot
iobj = m.Intermediate(m.integral(objective))

p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

# Equations
m.Equation(k.dt() == i*k - c)
# Objective Function
m.Maximize(final*m.integral(objective))

m.options.IMODE = 6
m.solve()

plt.figure(1)
plt.subplot(3,1,1)
plt.plot(m.time,c.value,'k:',linewidth=2,label=r'$C_{gekko}$')
C_sol = r*10*np.exp((i-r)*m.time)/(1-np.exp(-r*10))
plt.plot(m.time,C_sol,'r--',linewidth=2,label=r'$C_{exact}$')
plt.ylabel('Value'); plt.legend(loc='best')

plt.subplot(3,1,2)
plt.plot(m.time,k.value,'b-',linewidth=2,label=r'$K$')
plt.legend(loc='best')

plt.subplot(3,1,3)
plt.plot(m.time,objective.value,'g:',linewidth=2,label=r'$obj$')
plt.plot(m.time,iobj.value,'k',linewidth=2,label=r'$int obj$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()

Is there additional information that this problem is missing?

Edit: Added additional constraint k>0.

Added additional constraint as suggested in the comment. There is a small difference at the end from the exact solution because the last c value does not appear to influence the solution.

Solution with constraint

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

m = GEKKO(remote=True)
nt = 101; m.time = np.linspace(0,10,nt)
r = 0.33
i = 0.5

# Variables
c = m.MV(4,lb=0.001,ub=100); c.STATUS=1; c.DCOST=1e-6
m.free_initial(c)
k = m.Var(value=10,lb=0)
objective = m.Var(0)
t = m.Param(m.time)
m.Equation(objective==m.exp(-r*t)*m.log(c))

# just to include on the plot
iobj = m.Intermediate(m.integral(objective))

p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

# Equations
m.Equation(k.dt() == i*k - c)
# Objective Function
m.Maximize(final*m.integral(objective))

m.options.IMODE = 6
m.options.NODES = 3
m.solve()

plt.figure(1)
plt.subplot(3,1,1)
plt.plot(m.time,c.value,'k:',linewidth=2,label=r'$C_{gekko}$')
C_sol = r*10*np.exp((i-r)*m.time)/(1-np.exp(-r*10))
plt.plot(m.time,C_sol,'r--',linewidth=2,label=r'$C_{exact}$')
plt.ylabel('Value'); plt.legend(loc='best')

plt.subplot(3,1,2)
plt.plot(m.time,k.value,'b-',linewidth=2,label=r'$K$')
plt.legend(loc='best')

plt.subplot(3,1,3)
plt.plot(m.time,objective.value,'g:',linewidth=2,label=r'$obj$')
plt.plot(m.time,iobj.value,'k',linewidth=2,label=r'$int obj$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.show()

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome to Vigges Developer Community for programmer and developer-Open, Learning and Share
...