GEKKO - 使用自定义目标函数的参数估计 - 错误代码 -13 [英] GEKKO - parameter estimation with custom objective function - error code -13

查看:44
本文介绍了GEKKO - 使用自定义目标函数的参数估计 - 错误代码 -13的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

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

I have been successful in performing a steady-state parameter estimation employing the same techniques presented in the Gekko tutorials (linear and non-linear regression). Below is the code:

# -*- 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,Dy,Ho,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,0.52613,7.89266
0.46769,2.92548,0.53659,8.00228
0.45692,2.92548,0.5351,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,Dy,Ho,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,2.1979,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).我希望有人能指出我在哪里犯了错误.我认为目标函数是这样的,应该找到类似的结果.

However, when I try and implement a specific objective function I encounter an error (-13). I have tried setting a lower bound for the Var(s) however I then receive a different error code of (-1). I hope that that somebody can identify where I have made an error. I am of the view that the objective function is such that a similar result should be found.

# -*- 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,Dy,Ho,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,0.52613,7.89266
0.46769,2.92548,0.53659,8.00228
0.45692,2.92548,0.5351,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,Dy,Ho,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,2.1979,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_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 描述:

Error code -13 description for solver IPOPT is reported in the output:

EXIT: Invalid number in NLP function or derivative detected.

 An error occured.
 The error code is  -13

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

This means that there is likely divide by zero somewhere in the model. The other error -1 is when an infeasible solution is detected. The additional constraints prevent the solver from finding a solution because the equations and variable constraints can't be satisfied.

要找到成功的解决方案,需要做两件事:

Two things are needed to find a successful solution:

  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. 切换到 APOT 求解器:

m.options.SOLVER=1

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

There is now a successful solution.

# -*- 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,Dy,Ho,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,0.52613,7.89266
0.46769,2.92548,0.53659,8.00228
0.45692,2.92548,0.5351,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,Dy,Ho,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,2.1979,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

这篇关于GEKKO - 使用自定义目标函数的参数估计 - 错误代码 -13的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆