更改参数、目标函数和初始条件时使用 Gekko 优化套件出现“未找到解决方案”错误

【问题标题】："Solution Not Found" error using Gekko optimization suite when changing parameters, objective function, and initial condition更改参数、目标函数和初始条件时使用 Gekko 优化套件出现“未找到解决方案”错误
【发布时间】：2021-07-02 18:07:15
【问题描述】：

我正在尝试使用 Gekko 套件解决具有各种初始条件的多个优化问题。分配初始条件，使用 Gekko 运行优化，并收集每个解决方案。当我尝试更改参数、目标函数或初始条件时，Gekko 经常给我“未找到解决方案错误：第 2130 行，解决引发异常（apm_error）。”我在下面介绍一些案例，希望得到解决我遇到的这个问题的建议。我之前发布了一个类似的问题，但我希望这个问题更简洁明了。谢谢。

案例 1. 运行良好，没有错误。

from gekko import GEKKO
import pandas as pd
import numpy as np

dat = {'A0': [23221, 2198, 4296, 104906, 691], 'h': [0.04, 0.25, 0.07, 0, 12.58],'emax': [23221, 2198, 4296, 104906, 691] }
dftemp = pd.DataFrame(data=dat)
na=len(dftemp)

# time points
n=51
year=50

# constants
Pa0 = 3.061 
Pe0 = 10.603 
C0 = 100 
r = 0.05 
k=50 
shift=10000000 # to make positive inside log function
ll=0.15

  
for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=False)

    m.time = np.linspace(0,year,n)
    t=m.time
    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    
    
    A = m.SV(value=A0, lb=0, ub=A0)
    E = m.SV(value=0, lb=0, ub=A0)

    u = m.MV(value=0, lb=-emax, ub=emax) 
    u.STATUS = 1

    t = m.Param(value=m.time)
    C = m.Var(value=C0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(C.dt()==-C/k)
    m.Equation(d==m.exp(-t*r))
    # Objective (Utility)
    J = m.Var(value=0)
    
    # Final objective
    Jf = m.FV()
    Jf.STATUS = 1
    m.Connection(Jf,J,pos2='end')
    m.Equation(J.dt() == m.log((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift)*d)
    
    # maximize U
    m.Maximize(Jf)
    
    # options
    m.options.IMODE = 6  # optimal control
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 3 # solver (IPOPT)
    
    # solve optimization problem
    m.solve()
    
    # print profit
    print('Optimal Profit: ' + str(Jf.value[0]))

案例2.将“时间点”的数量从n=51更改为n=501发生错误

from gekko import GEKKO
import pandas as pd
import numpy as np

dat = {'A0': [23221, 2198, 4296, 104906, 691], 'h': [0.04, 0.25, 0.07, 0, 12.58],'emax': [23221, 2198, 4296, 104906, 691] }
dftemp = pd.DataFrame(data=dat)
na=len(dftemp)

# time points
n=501 
year=50

# constants
Pa0 = 3.061 
Pe0 = 10.603 
C0 = 100 
r = 0.05 
k=50 
shift=10000000 # to make positive inside log function
ll=0.15

   
for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=False)

    m.time = np.linspace(0,year,n)
    t=m.time
    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    
    
    A = m.SV(value=A0, lb=0, ub=A0)
    E = m.SV(value=0, lb=0, ub=A0)

    u = m.MV(value=0, lb=-emax, ub=emax) 
    u.STATUS = 1

    t = m.Param(value=m.time)
    C = m.Var(value=C0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(C.dt()==-C/k)
    m.Equation(d==m.exp(-t*r))
    # Objective (Utility)
    J = m.Var(value=0)
    
    # Final objective
    Jf = m.FV()
    Jf.STATUS = 1
    m.Connection(Jf,J,pos2='end')
    m.Equation(J.dt() == m.log((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift)*d)
    
    # maximize U
    m.Maximize(Jf)
    
    # options
    m.options.IMODE = 6  # optimal control
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 3 # solver (IPOPT)
    
    # solve optimization problem
    m.solve()
    
    # print profit
    print('Optimal Profit: ' + str(Jf.value[0]))

案例 3. 将目标函数从 m.log 更改为简单的线性求和。效果很好。

from gekko import GEKKO
import pandas as pd
import numpy as np

dat = {'A0': [23221, 2198, 4296, 104906, 691], 'h': [0.04, 0.25, 0.07, 0, 12.58],'emax': [23221, 2198, 4296, 104906, 691] }
dftemp = pd.DataFrame(data=dat)
na=len(dftemp)

# time points
n=51
year=50

# constants
Pa0 = 3.061 
Pe0 = 10.603 
C0 = 100 
r = 0.05 
k=50 
shift=10000000 # to make positive inside log function
ll=0.15
    
for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=False)

    m.time = np.linspace(0,year,n)
    t=m.time
    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    
    
    A = m.SV(value=A0, lb=0, ub=A0)
    E = m.SV(value=0, lb=0, ub=A0)

    u = m.MV(value=0, lb=-emax, ub=emax) 
    u.STATUS = 1

    t = m.Param(value=m.time)
    C = m.Var(value=C0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(C.dt()==-C/k)
    m.Equation(d==m.exp(-t*r))
    # Objective (Utility)
    J = m.Var(value=0)
    
    # Final objective
    Jf = m.FV()
    Jf.STATUS = 1
    m.Connection(Jf,J,pos2='end')
    m.Equation(J.dt() == ((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift)*d)
    
    # maximize U
    m.Maximize(Jf)
    
    # options
    m.options.IMODE = 6  # optimal control
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 3 # solver (IPOPT)
    
    # solve optimization problem
    m.solve()
    
    # print profit
    print('Optimal Profit: ' + str(Jf.value[0]))

案例 4. 将目标函数从 m.log 更改为简单求和，并从目标函数中删除“变量”d。发生错误

from gekko import GEKKO
import pandas as pd
import numpy as np

dat = {'A0': [23221, 2198, 4296, 104906, 691], 'h': [0.04, 0.25, 0.07, 0, 12.58],'emax': [23221, 2198, 4296, 104906, 691] }
dftemp = pd.DataFrame(data=dat)
na=len(dftemp)

# time points
n=51
year=50

# constants
Pa0 = 3.061 
Pe0 = 10.603 
C0 = 100 
r = 0.05 
k=50 
shift=10000000 # to make positive inside log function
ll=0.15

    
for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=False)

    m.time = np.linspace(0,year,n)
    t=m.time
    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    
    
    A = m.SV(value=A0, lb=0, ub=A0)
    E = m.SV(value=0, lb=0, ub=A0)

    u = m.MV(value=0, lb=-emax, ub=emax) 
    u.STATUS = 1

    t = m.Param(value=m.time)
    C = m.Var(value=C0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(C.dt()==-C/k)
    m.Equation(d==m.exp(-t*r))
    # Objective (Utility)
    J = m.Var(value=0)
    
    # Final objective
    Jf = m.FV()
    Jf.STATUS = 1
    m.Connection(Jf,J,pos2='end')
    m.Equation(J.dt() == ((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift))
    
    # maximize U
    m.Maximize(Jf)
    
    # options
    m.options.IMODE = 6  # optimal control
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 3 # solver (IPOPT)
    
    # solve optimization problem
    m.solve()
    
    # print profit
    print('Optimal Profit: ' + str(Jf.value[0]))

案例5.将目标函数从m.log改为简单的线性求和，将shift改为0，出现错误

from gekko import GEKKO
import pandas as pd
import numpy as np

dat = {'A0': [23221, 2198, 4296, 104906, 691], 'h': [0.04, 0.25, 0.07, 0, 12.58],'emax': [23221, 2198, 4296, 104906, 691] }
dftemp = pd.DataFrame(data=dat)
na=len(dftemp)

# time points
n=51
year=50

# constants
Pa0 = 3.061 
Pe0 = 10.603 
C0 = 100 
r = 0.05 
k=50 
shift=0 
ll=0.15
    
for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=False)

    m.time = np.linspace(0,year,n)
    t=m.time
    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    
    
    A = m.SV(value=A0, lb=0, ub=A0)
    E = m.SV(value=0, lb=0, ub=A0)

    u = m.MV(value=0, lb=-emax, ub=emax) 
    u.STATUS = 1

    t = m.Param(value=m.time)
    C = m.Var(value=C0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(C.dt()==-C/k)
    m.Equation(d==m.exp(-t*r))
    # Objective (Utility)
    J = m.Var(value=0)
    
    # Final objective
    Jf = m.FV()
    Jf.STATUS = 1
    m.Connection(Jf,J,pos2='end')
    m.Equation(J.dt() == ((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift)*d)
    
    # maximize U
    m.Maximize(Jf)
    
    # options
    m.options.IMODE = 6  # optimal control
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 3 # solver (IPOPT)
    
    # solve optimization problem
    m.solve()
    
    # print profit
    print('Optimal Profit: ' + str(Jf.value[0]))

【问题讨论】：

    标签： python optimization gekko


    【解决方案1】：

    对于其中许多情况，进行一些更改可能会有所帮助，例如：

    • 重新表述以避免m.log() 使用否定。中间变量也有帮助。
    # Objective (Utility)
J = m.Var(value=0)
rhs = m.Intermediate((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift)
m.Equation(m.exp(J.dt()/d)==rhs)
    • 在优化前初始化。 APOPT 求解器在初始化方面做得更好，IPOPT 在来自APOPT 的初始化解决方案上做得更好。初始化策略在article 中进行了讨论：Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E.，动态系统优化的初始化策略，计算机和化学工程，2015 年，卷。 78，第 39-50 页，DOI：10.1016/j.compchemeng.2015.04.016。
    # solve optimization problem
m.options.NODES = 3  # collocation nodes
m.options.SOLVER = 1 # solver (APOPT)
m.options.IMODE=4
m.solve()

# solve optimization problem
m.options.TIME_SHIFT=0
m.options.SOLVER=3 # solver (IPOPT)
m.options.IMODE=6
m.solve()

    案例 2：n=501 的成功解决方案

    from gekko import GEKKO
import pandas as pd
import numpy as np

dat = {'A0': [23221, 2198, 4296, 104906, 691], \
       'h': [0.04, 0.25, 0.07, 0, 12.58],\
       'emax': [23221, 2198, 4296, 104906, 691] }
dftemp = pd.DataFrame(data=dat)
na=len(dftemp)

# time points
n=501 
year=50

# constants
Pa0 = 3.061 
Pe0 = 10.603 
C0 = 100 
r = 0.05 
k=50 
shift=10000000 # to make positive inside log function
ll=0.15

for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=True)

    m.time = np.linspace(0,year,n)
    t=m.time
    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    
    
    A = m.SV(value=A0, lb=0, ub=A0)
    E = m.SV(value=0, lb=0, ub=A0)

    u = m.MV(value=0, lb=-emax, ub=emax) 
    u.STATUS = 1

    t = m.Param(value=m.time)
    C = m.Var(value=C0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(C.dt()==-C/k)
    m.Equation(d==m.exp(-t*r))
    # Objective (Utility)
    J = m.Var(value=0)
    rhs = m.Intermediate((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift)
    m.Equation(m.exp(J.dt()/d)==rhs)
        
    # Final objective
    final = np.zeros_like(m.time)
    final[-1] = 1
    final = m.Param(final)
    
    # maximize U
    m.Maximize(final*J)
        
    # solve optimization problem
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 1 # solver (APOPT)
    m.options.IMODE=4
    print('\n\nInitializing with APOPT')
    m.solve(disp=False)

    # solve optimization problem
    m.options.TIME_SHIFT=0
    m.options.SOLVER=3 # solver (IPOPT)
    m.options.IMODE=6
    m.solve()
        
    # print profit
    print('Optimal Profit: ' + str(J.value[-1]))

    请在其他情况下也尝试这种方法。通常使用 m.sum() 而不是 sum 更好，因为 Gekko 将求和构造为内置对象而不是长字符串的方式。

    编辑：更新的解决方案这是具有正确初始条件的更新解决方案。 IMODE=4 的初始条件不会转移到 IMODE=6。另一种方法是将IMODE=6 与COLDSTART=1 等效使用。

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

dat = {'A0': [23221, 2198, 4296, 104906, 691], \
       'h': [0.04, 0.25, 0.07, 0, 12.58],\
       'emax': [23221, 2198, 4296, 104906, 691] }
dftemp = pd.DataFrame(data=dat)
na=len(dftemp)

# time points
n=101 
year=50

# constants
Pa0 = 3.061 
Pe0 = 10.603 
C0 = 100 
r = 0.05 
k=50 
shift=10000000 # to make positive inside log function
ll=0.15

plt.figure(figsize=(10,5))
for i in range(0,na):
    # create GEKKO model
    m = GEKKO(remote=True)

    m.time = np.linspace(0,year,n)
    t=m.time
    
    A0=dftemp.loc[i][0]
    h=dftemp.loc[i][1]
    emax=dftemp.loc[i][2]    

    print('A (initial): ' + str(A0))
    A = m.Var(value=A0, lb=0, ub=A0)
    E = m.Var(value=0, lb=0, ub=A0)

    u = m.MV(value=0, lb=-emax, ub=emax) 
    u.STATUS = 1

    t = m.Param(value=m.time)
    C = m.Var(value=C0)
    d = m.Var(value=1)
    l = m.Param(value=ll)
    
    # Equation
    m.Equation(A.dt()==-u)
    m.Equation(E.dt()==u)
    m.Equation(C.dt()==-C/k)
    m.Equation(d==m.exp(-t*r))
    # Objective (Utility)
    J = m.Var(value=0)
    rhs = m.Intermediate((A+E*(1-l))*h*Pa0-C*u+E*Pe0+shift)
    m.Equation(m.exp(J.dt()/d)==rhs)
        
    # Final objective
    final = np.zeros_like(m.time)
    final[-1] = 1
    final = m.Param(final)
    
    # maximize U
    m.Maximize(final*J)
        
    # solve optimization problem
    m.options.NODES = 3  # collocation nodes
    m.options.SOLVER = 1 # solver (APOPT)
    m.options.IMODE=6
    m.options.COLDSTART=1
    print('\n\nInitializing with APOPT')
    m.solve(disp=False)

    # solve optimization problem
    m.options.COLDSTART=0
    m.options.TIME_SHIFT=0
    m.options.SOLVER=3 # solver (IPOPT)
    m.options.IMODE=6
    m.solve()
        
    # print profit
    print('Optimal Profit: ' + str(J.value[-1]))

    plt.plot(m.time,A.value,label='A case '+str(i))
    plt.plot(m.time,E.value,label='E case '+str(i))
    plt.plot(m.time,u.value,label='u case '+str(i))
    

plt.legend()
plt.show()

    【讨论】：

    • 感谢您的回答@JohnHedengren。我想问一下解决方案。恐怕答案中的解决方案与我的预期有所不同。答案代码似乎使用了固定的初始条件； 1 用于初始A 和E，这与我想要使用的不同（来自 dat 数据帧的A0s）。你能告诉我这里缺少什么吗？
    • 我解决了这个问题。请查看编辑。
    • 谢谢@JohnHedengren。 A 和 E 的问题已解决。我认为此时，u 似乎以某种方式增加了额外的限制。我期待解决方案在特定时间出现u 的峰值，但u 在整个时间内均匀分布，其上限远小于我设置的（或使A 为零的大小）在最后时间）。
    • 我可以再问一个问题吗？如果我将u 的上限从emax 替换为A.value，我是否正确理解这使得时间(t) 的上限变为上一次时间(t-1) 的A.value？我想要做的是将u 的边界限制为状态变量的最大值，每次都是A 和E。
    • 您可以设置一个新的不等式约束，例如m.Equation(u<=A) 和m.Equation(u<=E)。如果上限是A 和E 的最大值，则可以使用m.Equation(u<=m.max3(A,E))。定义u.UPPER = max(A.value) 只使用初始化时A.value 的当前值，不会给出想要的结果。
    猜你喜欢
    • 1970-01-01
    • 2022-01-14
    • 1970-01-01
    • 2022-10-14
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2011-09-04
    • 2016-08-19
    相关资源
    最近更新 更多
    热门标签
    Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode