title: GEKKO - parameter estimation with custom objective function - error code -13
body: <p>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:</p>
  <pre class="lang-py prettyprint-override"><code># -*- coding: utf-8 -*-
  &quot;&quot;&quot;
  Spyder Editor
  
  This is a temporary script file.
  &quot;&quot;&quot;
  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 = [&quot;Tb&quot;, &quot;Dy&quot;, &quot;Ho&quot;, &quot;Y&quot;]
  
  # import dataframe
  
  x_i_a_Lns_v = (
      pd.read_csv(
          StringIO(
              &quot;&quot;&quot;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&quot;&quot;&quot;
          )
      )
      / 1000
  )
  
  x_i_a_H_v = pd.read_csv(
      StringIO(
          &quot;&quot;&quot;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&quot;&quot;&quot;
      )
  )
  
  x_i_o_HA2_v = pd.read_csv(
      StringIO(
          &quot;&quot;&quot;(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&quot;&quot;&quot;
      )
  )
  
  x_f_a_Lns_v = (
      pd.read_csv(
          StringIO(
              &quot;&quot;&quot;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&quot;&quot;&quot;
          )
      )
      / 1000
  )
  
  # Issue with NaN (missing value)
  idx = x_f_a_Lns_v[x_f_a_Lns_v[&quot;Ho&quot;] &gt;= 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, &quot;H&quot;].values)
  x_i_o_HA2 = m.Param(value=x_i_o_HA2_v.loc[idx, &quot;(HA)2&quot;].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&quot;K_eq_{v}: {K_eqs[i][0]}&quot;)
  
  </code></pre>
  <p>Output:</p>
  <pre><code>K_eq_Tb: 5.7327051238
  K_eq_Dy: 15.581534791
  K_eq_Ho: 27.414244408
  K_eq_Y: 46.925190325
  </code></pre>
  <p>This is as expected.</p>
  <p>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.</p>
  <pre class="lang-py prettyprint-override"><code># -*- coding: utf-8 -*-
  &quot;&quot;&quot;
  Spyder Editor
  
  This is a temporary script file.
  &quot;&quot;&quot;
  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 = [&quot;Tb&quot;, &quot;Dy&quot;, &quot;Ho&quot;, &quot;Y&quot;]
  
  # import dataframes
  x_i_a_Lns_v = (
      pd.read_csv(
          StringIO(
              &quot;&quot;&quot;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&quot;&quot;&quot;
          )
      )
      / 1000
  )
  
  x_i_a_H_v = pd.read_csv(
      StringIO(
          &quot;&quot;&quot;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&quot;&quot;&quot;
      )
  )
  
  x_i_o_HA2_v = pd.read_csv(
      StringIO(
          &quot;&quot;&quot;(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&quot;&quot;&quot;
      )
  )
  
  x_f_a_Lns_v = (
      pd.read_csv(
          StringIO(
              &quot;&quot;&quot;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&quot;&quot;&quot;
          )
      )
      / 1000
  )
  
  # Issue with NaN (missing value)
  idx = x_f_a_Lns_v[x_f_a_Lns_v[&quot;Ho&quot;] &gt;= 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, &quot;H&quot;].values)
  x_i_o_HA2 = m.Param(value=x_i_o_HA2_v.loc[idx, &quot;(HA)2&quot;].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&quot;K_eq_{v}: {K_eqs[i][0]}&quot;)
  
  </code></pre>