微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦！

GEKKO - 使用自定义目标函数的参数估计 - 错误代码 -13

分类：编程问答

如何解决GEKKO - 使用自定义目标函数的参数估计 - 错误代码 -13

我使用 Gekko 教程中介绍的相同技术（线性和非线性回归）成功地执行了稳态参数估计。代码如下：

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
from io import StringIO

from gekko import GEKKO

import numpy as np
import pandas as pd


import matplotlib.pyplot as plt

# Tb,Dy,Ho,Y,H,(HA)2
# measurements of the initial and final aqueous concentrations

aqueous_species = ["Tb","Dy","Ho","Y"]

# import dataframe

x_i_a_Lns_v = (
    pd.read_csv(
        StringIO(
            """Tb,Y
0.19538,1.22628,0.2242,3.39823
0.28462,1.83411,0.32435,4.90551
0.34769,2.1979,0.39609,5.89209
0.37692,2.39495,0.43794,6.57722
0.41231,2.62232,0.47382,7.09791
0.43538,2.78906,0.5052,7.45418
0.44462,2.88,0.51866,7.89266
0.46,2.92548,0.52912,7.83785
0.46615,2.94064,0.5351,7.97488
0.45846,2.91032,0.52613,7.94747
0.46,2.86485,7.89266
0.46769,0.53659,8.00228
0.45692,8.02969
0.47385,2.98611,0.55005,8.13931
0.47385,2.97095,0.54108,8.1119"""
        )
    )
    / 1000
)

x_i_a_H_v = pd.read_csv(
    StringIO(
        """H
10.01809
7.28346
5.16795
3.62036
2.50218
1.7173
1.17411
0.80491
0.54616
0.37078
0.25406
0.16932
0.11262
0.07574
0.04455"""
    )
)

x_i_o_HA2_v = pd.read_csv(
    StringIO(
        """(HA)2
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746"""
    )
)

x_f_a_Lns_v = (
    pd.read_csv(
        StringIO(
            """Tb,Y
0.19231,1.23234,0.21822,3.34342
0.26923,1.74316,0.31239,4.60405
0.33692,2.13727,0.38862,5.72766
0.37385,2.33432,0.41253,6.05652
0.40923,2.50106,0.43196,5.97431
0.38769,0.3363,3.91892
0.33692,1.53095,0.19431,1.91287
0.24,0.82307,0.0852,0.77282
0.12,0.33347,0.03139,0.27679
0.05846,0.14552,0.01495,0.1151
0.02615,0.06669,0.00598,0.05207
0.01231,0.03032,0.02466,0.00548
0.00615,0.01364,0.01096,0.00548"""
        )
    )
    / 1000
)

# Issue with NaN (missing value)
idx = x_f_a_Lns_v[x_f_a_Lns_v["Ho"] >= 0].index
# idx = x_f_a_Lns_v.index

#### Solution
# create model to be solved locally
m = GEKKO(remote=False)

x_i_a_Lns = [m.Param(value=x_i_a_Lns_v.loc[idx,v].values) for v in aqueous_species]

x_i_a_H = m.Param(value=x_i_a_H_v.loc[idx,"H"].values)
x_i_o_HA2 = m.Param(value=x_i_o_HA2_v.loc[idx,"(HA)2"].values)

x_f_a_Lns = [m.CV(value=x_f_a_Lns_v.loc[idx,v].values) for v in aqueous_species]

for i in range(len(aqueous_species)):
    x_f_a_Lns[i].FSTATUS = 1

# equilibrium constants
K_eqs = [m.FV(value=1) for i in range(len(aqueous_species))]

for i in range(len(x_f_a_Lns_v.columns)):
    K_eqs[i].STATUS = 1


# equilibrium model
x_f_a_H = m.Intermediate(
    x_i_a_H
    + 3 * m.sum([(x_i_a_Lns[i] - x_f_a_Lns[i]) for i in range(len(aqueous_species))])
)
x_f_o_HA2 = m.Intermediate(
    x_i_o_HA2
    - 3 * m.sum([(x_i_a_Lns[i] - x_f_a_Lns[i]) for i in range(len(aqueous_species))])
)

m.Equations(
    [
        K_eqs[i]
        == ((x_i_a_Lns[i] - x_f_a_Lns[i]) * x_f_a_H ** 3)
        / (x_f_a_Lns[i] * x_f_o_HA2 ** 3)
        for i in range(len(aqueous_species))
    ]
)


# regression
m.options.IMODE = 2

m.options.EV_TYPE = 2
m.solve(disp=True)

for i,v in enumerate(aqueous_species):
    print(f"K_eq_{v}: {K_eqs[i][0]}")

输出：

K_eq_Tb: 5.7327051238
K_eq_Dy: 15.581534791
K_eq_Ho: 27.414244408
K_eq_Y: 46.925190325

这符合预期。

然而，当我尝试实现一个特定的目标函数时，我遇到了一个错误 (-13)。我曾尝试为 Var(s) 设置下限，但随后我收到了不同的错误代码 (-1)。我希望有人能指出我在哪里犯了错误。我认为目标函数是这样的，应该找到类似的结果。

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
from io import StringIO

from gekko import GEKKO

import numpy as np
import pandas as pd


import matplotlib.pyplot as plt

# Tb,"Y"]

# import dataframes
x_i_a_Lns_v = (
    pd.read_csv(
        StringIO(
            """Tb,"(HA)2"].values)


x_f_a_Lns_m = [m.Param(value=x_f_a_Lns_v.loc[idx,v].values) for v in aqueous_species]
x_f_a_Lns = [m.Var() for v in aqueous_species]
# x_f_a_Lns = [m.Var(lb=1e-12) for v in aqueous_species]

# equilibrium constants
K_eqs = [m.FV(value=1) for i in range(len(aqueous_species))]

for i in range(len(aqueous_species)):
    K_eqs[i].STATUS = 1


# equilibrium model
x_f_a_H = m.Intermediate(
    x_i_a_H
    + 3 * m.sum([(x_i_a_Lns[i] - x_f_a_Lns[i]) for i in range(len(aqueous_species))])
)
x_f_o_HA2 = m.Intermediate(
    x_i_o_HA2
    - 3 * m.sum([(x_i_a_Lns[i] - x_f_a_Lns[i]) for i in range(len(aqueous_species))])
)

m.Equations(
    [
        K_eqs[i]
        == ((x_i_a_Lns[i] - x_f_a_Lns[i]) * x_f_a_H ** 3)
        / (x_f_a_Lns[i] * x_f_o_HA2 ** 3)
        for i in range(len(aqueous_species))
    ]
)


m.Obj(
    m.sum(
        [
            ((x_f_a_Lns_m[i] - x_f_a_Lns[i]) / x_f_a_Lns_m[i]) ** 2
            for i in range(len(aqueous_species))
        ]
    )
)

# regression
m.options.IMODE = 2

# m.options.EV_TYPE = 2
m.solve(disp=True)

for i,v in enumerate(aqueous_species):
    print(f"K_eq_{v}: {K_eqs[i][0]}")

解决方法

输出中报告了求解器IPOPT的错误代码-13描述：

EXIT: Invalid number in NLP function or derivative detected.

 An error occured.
 The error code is  -13

这意味着模型中的某处可能被零除。另一个错误 -1 是当检测到不可行的解决方案时。由于无法满足方程和变量约束，额外的约束会阻止求解器找到解。

要找到成功的解决方案，需要做两件事：

  1. 重新排列方程以避免被零除。
m.Equations(
    [
        K_eqs[i]
        == ((x_i_a_Lns[i] - x_f_a_Lns[i]) * x_f_a_H ** 3)
        / (x_f_a_Lns[i] * x_f_o_HA2 ** 3)
        for i in range(len(aqueous_species))
    ]
)

重新排列为：

m.Equations(
    [
        K_eqs[i] * (x_f_a_Lns[i] * x_f_o_HA2 ** 3)
        == ((x_i_a_Lns[i] - x_f_a_Lns[i]) * x_f_a_H ** 3)
        for i in range(len(aqueous_species))
    ]
)
  1. 切换到 APPT 求解器：
m.options.SOLVER=1

现在有一个成功的解决方案。

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
from io import StringIO

from gekko import GEKKO

import numpy as np
import pandas as pd


import matplotlib.pyplot as plt

# Tb,Dy,Ho,Y,H,(HA)2
# measurements of the initial and final aqueous concentrations

aqueous_species = ["Tb","Dy","Ho","Y"]

# import dataframes
x_i_a_Lns_v = (
    pd.read_csv(
        StringIO(
            """Tb,Y
0.19538,1.22628,0.2242,3.39823
0.28462,1.83411,0.32435,4.90551
0.34769,2.1979,0.39609,5.89209
0.37692,2.39495,0.43794,6.57722
0.41231,2.62232,0.47382,7.09791
0.43538,2.78906,0.5052,7.45418
0.44462,2.88,0.51866,7.89266
0.46,2.92548,0.52912,7.83785
0.46615,2.94064,0.5351,7.97488
0.45846,2.91032,0.52613,7.94747
0.46,2.86485,7.89266
0.46769,0.53659,8.00228
0.45692,8.02969
0.47385,2.98611,0.55005,8.13931
0.47385,2.97095,0.54108,8.1119"""
        )
    )
    / 1000
)

x_i_a_H_v = pd.read_csv(
    StringIO(
        """H
10.01809
7.28346
5.16795
3.62036
2.50218
1.7173
1.17411
0.80491
0.54616
0.37078
0.25406
0.16932
0.11262
0.07574
0.04455"""
    )
)

x_i_o_HA2_v = pd.read_csv(
    StringIO(
        """(HA)2
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746
0.4746"""
    )
)

x_f_a_Lns_v = (
    pd.read_csv(
        StringIO(
            """Tb,Y
0.19231,1.23234,0.21822,3.34342
0.26923,1.74316,0.31239,4.60405
0.33692,2.13727,0.38862,5.72766
0.37385,2.33432,0.41253,6.05652
0.40923,2.50106,0.43196,5.97431
0.38769,0.3363,3.91892
0.33692,1.53095,0.19431,1.91287
0.24,0.82307,0.0852,0.77282
0.12,0.33347,0.03139,0.27679
0.05846,0.14552,0.01495,0.1151
0.02615,0.06669,0.00598,0.05207
0.01231,0.03032,0.02466,0.00548
0.00615,0.01364,0.01096,0.00548"""
        )
    )
    / 1000
)

# Issue with NaN (missing value)
idx = x_f_a_Lns_v[x_f_a_Lns_v["Ho"] >= 0].index
# idx = x_f_a_Lns_v.index

#### Solution
# create model to be solved locally
m = GEKKO(remote=False)

x_i_a_Lns = [m.Param(value=x_i_a_Lns_v.loc[idx,v].values) for v in aqueous_species]

x_i_a_H = m.Param(value=x_i_a_H_v.loc[idx,"H"].values)
x_i_o_HA2 = m.Param(value=x_i_o_HA2_v.loc[idx,"(HA)2"].values)


x_f_a_Lns_m = [m.Param(value=x_f_a_Lns_v.loc[idx,v].values) for v in aqueous_species]
x_f_a_Lns = [m.Var() for v in aqueous_species]
# x_f_a_Lns = [m.Var(lb=1e-12) for v in aqueous_species]

# equilibrium constants
K_eqs = [m.FV(value=1) for i in range(len(aqueous_species))]

for i in range(len(aqueous_species)):
    K_eqs[i].STATUS = 1


# equilibrium model
x_f_a_H = m.Intermediate(
    x_i_a_H
    + 3 * m.sum([(x_i_a_Lns[i] - x_f_a_Lns[i]) for i in range(len(aqueous_species))])
)
x_f_o_HA2 = m.Intermediate(
    x_i_o_HA2
    - 3 * m.sum([(x_i_a_Lns[i] - x_f_a_Lns[i]) for i in range(len(aqueous_species))])
)

m.Equations(
    [
        K_eqs[i] * (x_f_a_Lns[i] * x_f_o_HA2 ** 3)
        == ((x_i_a_Lns[i] - x_f_a_Lns[i]) * x_f_a_H ** 3)
        for i in range(len(aqueous_species))
    ]
)


m.Obj(
    m.sum(
        [
            ((x_f_a_Lns_m[i] - x_f_a_Lns[i]) / x_f_a_Lns_m[i]) ** 2
            for i in range(len(aqueous_species))
        ]
    )
)

# regression
m.options.IMODE = 2

m.options.SOLVER=1

# m.options.EV_TYPE = 2
m.solve(disp=True)

for i,v in enumerate(aqueous_species):
    print(f"K_eq_{v}: {K_eqs[i][0]}")

 Successful solution
 
 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.06899999999999999 sec
 Objective      :  0.32417750846391324
 Successful solution
 ---------------------------------------------------
 

K_eq_Tb: 5.489166594
K_eq_Dy: 15.245134386
K_eq_Ho: 30.529918885
K_eq_Y: 54.83667882

版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。

苹果市值2025年有望达4万亿美元
pythonJavaScriptjavaHTMLPHPreactjsC#AndroidCSSNode.jssqlrpython-3.xMysqLjQueryc++pandasFlutterangularIOSdjangolinuxswifttypescript路由器JSON路由器设置无线路由器h3c华三华三路由器设置华三路由器电脑软件教程arraysdocker软件图文教程Cvue.jslaravelspring-boot