2 votes

Kurt Forrester avatar

GEKKO - estimation des paramètres avec une fonction objective personnalisée - code d'erreur -13

Demandé el 18 de Avril, 2021
Quand la question a-t-elle été
37 affichage
Nombre de visites la question a
1 Réponses
Nombre de réponses aux questions
Résolu
Situation réelle de la question

J'ai réussi à effectuer une estimation des paramètres en régime permanent en utilisant les mêmes techniques que celles présentées dans les tutoriels de Gekko (régression linéaire et non linéaire). Voici le 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]}")

Sortie :

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

C'est ce que l'on attendait.

Cependant, lorsque j'essaie de mettre en œuvre une fonction objective spécifique, je rencontre une erreur (-13). J'ai essayé de définir une limite inférieure pour la (les) Var(s), mais je reçois alors un autre code d'erreur (-1). J'espère que quelqu'un pourra me dire où j'ai fait une erreur. Je pense que la fonction objectif est telle qu'un résultat similaire devrait être trouvé.

# -*- 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]}")
Demandé el 18 de Avril, 2021 par Kurt Forrester

Réponse

Trop de publicités?

1voto

John Hedengren avatar
John Hedengren Points 1

La description du code d'erreur -13 pour le solveur IPOPT est signalée dans la sortie :

EXIT: Invalid number in NLP function or derivative detected.

 An error occured.
 The error code is  -13

Cela signifie qu'il y a probablement une division par zéro quelque part dans le modèle. L'autre erreur -1 se produit lorsqu'une solution infaisable est détectée. Les contraintes supplémentaires empêchent le solveur de trouver une solution parce que les équations et les contraintes variables ne peuvent pas être satisfaites.

Deux éléments sont nécessaires pour trouver une solution efficace :

  1. Réarrangez l'équation pour éviter de diviser par zéro.

    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)) ] )

Réarranger en :

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. Passer au solveur APOPT :

    m.options.SOLVER=1

Il existe désormais une solution efficace.

# -*- 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
Répondu el 19 de Avril, 2021 par John Hedengren (1 Points )

Questions en vedette

Top Tags

Dans notre réseau

Prograide.com

Prograide est une communauté de développeurs qui cherche à élargir la connaissance de la programmation au-delà de l'anglais.
Pour cela nous avons les plus grands doutes résolus en français et vous pouvez aussi poser vos propres questions ou résoudre celles des autres.

Powered by:

Yandex
X