How to interpret shadow price array shape in Gekko

Question

I'm trying to calculate the shadow prices for the various constraints and place them in a dictionary following the form c[i]:sp[i] where c is the name of a specific constraint and sp is the numeric shadow price value. I am able to get "apm_lam.txt" to generate locally, but the shape of the of the array is 506, which is way more than the number of constraints that should have been added to the model (I'm counting like 60-something).

Is there an easy to tie-back the shadow prices in this output file to the actual named constraints in the Gekko model?

import numpy as np
import pandas as pd
from gekko import GEKKO
m = GEKKO(remote=False)
m.options.NODES = 3
m.options.IMODE = 3
m.options.MAX_ITER = 1000

lnuc_weeks = [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]

min_promo_price = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,3]

max_promo_price = [3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5,3.5, 3.5, 3.5, 3.5, 3.5, 3.5]

base_srp = [3.48, 3.48, 3.48, 3.48, 3.0799, 3.0799, 3.0799, 3.0799,3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799, 3.0799]

lnuc_min_promo_price = 1.99

lnuc_max_promo_price = 1.99

coeff_fedi = [0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589,0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589, 0.022589]

coeff_feao = [0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995, 0.02929995]

coeff_diso = [0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338, 0.05292338]

sumproduct_base = [0.20560305, 0.24735297, 0.24957423, 0.23155435, 0.23424058,0.2368096 , 0.27567109, 0.27820648, 0.2826393 , 0.28660598, 0.28583971, 0.30238505, 0.31726649, 0.31428312, 0.31073792, 0.29036779, 0.32679041, 0.32156337, 0.24633734]

neg_ln = [[0.14842000515],[0.14842000512],[0.14842000515],[0.14842000512],[-0.10407483058],[0.43676249024],[0.43676249019],[0.43676249024],[0.43676249019],[0.43676249024],[0.43676249019], [0.026284840258],[0.026284840291],[0.026284840258],[0.026284840291], [0.026185109811],[0.026284840258],[0.026284840291],[0.026284840258]]

neg_ln_ppi_coeff = [1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879,1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879, 1.22293879,1.22293879, 1.22293879, 1.22293879, 1.22293879]

base_volume = [124.38, 193.2, 578.72, 183.88, 197.42, 559.01, 67.68, 110.01,60.38, 177.11, 102.65, 66.02, 209.83, 81.22, 250.44, 206.44, 87.99, 298.95, 71.07]

week = pd.Series([13, 14, 17, 18, 19, 26, 28, 33, 34, 35, 39, 42, 45, 46, 47, 48, 50, 51, 52])


n = 19

cons = []
shadow = []

x1 = m.Array(m.Var,(n), integer=True) #LNUC weeks

i = 0
for xi in x1:
    xi.value = lnuc_weeks[i]
    xi.lower = 0
    xi.upper = lnuc_weeks[i]
    i += 1

x2 = m.Array(m.Var,(n)) #Blended SRP

i = 0
for xi in x2:
    xi.value = 5
    m.Equation(xi >= m.if3((x1[i]) - 0.5, min_promo_price[i], lnuc_min_promo_price))
    m.Equation(xi <= m.if3((x1[i]) - 0.5, max_promo_price[i], lnuc_max_promo_price))
    i += 1
    
x3 = m.Array(m.Var,(n), integer=True) #F&D
x4 = m.Array(m.Var,(n), integer=True) #FO
x5 = m.Array(m.Var,(n), integer=True) #DO
x6 = m.Array(m.Var,(n), integer=True) #TPR

#Default to F&D
i = 0
for xi in x3:
    xi.value = 1
    xi.lower = 0
    xi.upper = 1
    i += 1

i = 0
for xi in x4:
    xi.value = 0
    xi.lower = 0
    xi.upper = 1
    i += 1

i = 0
for xi in x5:
    xi.value = 0
    xi.lower = 0
    xi.upper = 1
    i += 1

i = 0
for xi in x6:
    xi.value = 0
    xi.lower = 0
    xi.upper = 1
    i += 1

x7 = m.Array(m.Var,(n), integer=True) #Max promos

i = 0
for xi in x7:
    xi.value = 1
    xi.lower = 0
    xi.upper = 1
    i += 1

x = [x1,x2,x3,x4,x5,x6,x7]

neg_ln=[m.Intermediate(-m.log(x[1][i]/base_srp[i])) for i in range(n)]

total_vol_fedi  =[m.Intermediate(coeff_fedi[0]+ sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]
total_vol_feao  =[m.Intermediate(coeff_feao[0]+ sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]
total_vol_diso  =[m.Intermediate(coeff_diso[0]+ sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]
total_vol_tpro  =[m.Intermediate(sumproduct_base[i] + (neg_ln[i]*neg_ln_ppi_coeff[0])) for i in range(n)]

simu_total_volume = [m.Intermediate((
(m.max2(0,base_volume[i]*(m.exp(total_vol_fedi[i])-1)) * x[2][i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_feao[i])-1)) * x[3][i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_diso[i])-1)) * x[4][i] +
m.max2(0,base_volume[i]*(m.exp(total_vol_tpro[i])-1)) * x[5][i]) + base_volume[i]) * x[6][i]) for i in range(n)]


[m.Equation(x3[i] + x4[i] + x5[i] + x6[i] == 1) for i in range(i)]
    
    

#Limit max promos
m.Equation(sum(x7)<=10)

#Enforce spacing and duration
s=1
for s2 in range(1, s+1):
    for i in range(0, n-s2):
        f = week[week == week[i] + s2].index
        if len(f) > 0:
            m.Equation(x7[i] + x7[f[0]]<=1)

m.Maximize(m.sum(simu_total_volume))

m.options.SOLVER=3
m.options.DIAGLEVEL = 2
m.solve(disp = True, debug=True)

df = pd.concat([pd.Series(week), pd.Series([i[0] for i in x7]), pd.Series([i[0] for i in simu_total_volume])], axis=1)
df.columns = ['week', 'x7', 'total_volume']
df[df['x7']>0]

m.open_folder()
lam = np.loadtxt(m.path+'/apm_lam.txt')
print(lam.shape)

Answer 1

Additional information on the constraints is given in the rto_4_eqn.txt file in the run directory when remote=False. Add this code after m.solve() to view a report on the equations:

fid = open(m.path+'/rto_4_eqn.txt','r')
print(fid.read())
fid.close()

For different modes of operation, the first 3 characters of the file change.

IMODE	Mode	Equation Summary Filename
1	Steady-State Simulation	ss_4_eqn.txt
2	Model Parameter Update	mpu_4_eqn.txt
3	Real-Time Optimization	rto_4_eqn.txt
4	Dynamic Simulation	sim_4_eqn.txt
5	Dynamic Estimation	est_4_eqn.txt
6	Dynamic Control	ctl_4_eqn.txt
7	Sequential Simulation	sqs_4_eqn.txt
8	Sequential Estimation	est_4_eqn.txt
9	Sequential Control	ctl_4_eqn.txt

Look for the section on ACTIVE EQUATIONS with 506 equations.

************************************************
************* ACTIVE EQUATIONS *****************
************************************************
Number           ID  Node    Horizon  Unscaled Res  Scaled Res   Scaling     Name
         1          1    1          1   1.0000E-08   1.0000E-08   1.0000E+00  ss.max2_1.mpcc[1]: x[2] - x[1] = s[1] - s[2]
         2          2    1          1  -1.0000E-08  -1.0000E-08   1.0000E+00  ss.max2_1.output_eqn: y = x[2] + s[2]
         3          4    1          1   1.0000E-08   1.0000E-08   1.0000E+00  ss.max2_2.mpcc[1]: x[2] - x[1] = s[1] - s[2]
         4          5    1          1  -1.0000E-08  -1.0000E-08   1.0000E+00  ss.max2_2.output_eqn: y = x[2] + s[2]
         5          7    1          1   1.0000E-08   1.0000E-08   1.0000E+00  ss.max2_3.mpcc[1]: x[2] - x[1] = s[1] - s[2]
         6          8    1          1  -1.0000E-08  -1.0000E-08   1.0000E+00  ss.max2_3.output_eqn: y = x[2] + s[2]
         7         10    1          1   1.0000E-08   1.0000E-08   1.0000E+00  ss.max2_4.mpcc[1]: x[2] - x[1] = s[1] - s[2]
         8         11    1          1  -1.0000E-08  -1.0000E-08   1.0000E+00  ss.max2_4.output_eqn: y = x[2] + s[2]
         9         13    1          1   1.0000E-08   1.0000E-08   1.0000E+00  ss.max2_5.mpcc[1]: x[2] - x[1] = s[1] - s[2]
        10         14    1          1  -1.0000E-08  -1.0000E-08   1.0000E+00  ss.max2_5.output_eqn: y = x[2] + s[2]
...
       484        560    1          1  -1.1718E-12  -1.1718E-12   1.0000E+00  ss.Eqn(331): 0 = 1-((int_v204+int_v205))-slk_331
       485        561    1          1   0.0000E+00   0.0000E+00   1.0000E+00  ss.Eqn(332): 0 = 1-((int_v205+int_v206))-slk_332
       486        562    1          1   0.0000E+00   0.0000E+00   1.0000E+00  ss.Eqn(333): 0 = 1-((int_v207+int_v208))-slk_333
       487        563    1          1   0.0000E+00   0.0000E+00   1.0000E+00  ss.Eqn(334): 0 = 1-((int_v208+int_v209))-slk_334
       488        564    1          1  -1.8406E-06  -1.8406E-06   1.0000E+00  ss.Eqn(335): 0 = v438-(i209)
       489        565    1          1   2.2181E-06   2.2181E-06   1.0000E+00  ss.Eqn(336): 0 = v439-(i210)
       490        566    1          1   6.5732E-06   6.5732E-06   1.0000E+00  ss.Eqn(337): 0 = v440-(i211)
       491        567    1          1  -5.0622E-09  -5.0622E-09   1.0000E+00  ss.Eqn(338): 0 = v441-(i212)
       492        568    1          1   2.3812E-06   2.3812E-06   1.0000E+00  ss.Eqn(339): 0 = v442-(i213)
       493        569    1          1   6.6581E-06   6.6581E-06   1.0000E+00  ss.Eqn(340): 0 = v443-(i214)
       494        570    1          1  -1.5305E-06  -1.5305E-06   1.0000E+00  ss.Eqn(341): 0 = v444-(i215)
       495        571    1          1   1.3645E-06   1.3645E-06   1.0000E+00  ss.Eqn(342): 0 = v445-(i216)
       496        572    1          1  -1.3760E-06  -1.3760E-06   1.0000E+00  ss.Eqn(343): 0 = v446-(i217)
       497        573    1          1   2.1610E-06   2.1610E-06   1.0000E+00  ss.Eqn(344): 0 = v447-(i218)
       498        574    1          1  -1.1720E-06  -1.1720E-06   1.0000E+00  ss.Eqn(345): 0 = v448-(i219)
       499        575    1          1  -9.3201E-07  -9.3201E-07   1.0000E+00  ss.Eqn(346): 0 = v449-(i220)
       500        576    1          1   2.3803E-06   2.3803E-06   1.0000E+00  ss.Eqn(347): 0 = v450-(i221)
       501        577    1          1  -1.1581E-06  -1.1581E-06   1.0000E+00  ss.Eqn(348): 0 = v451-(i222)
       502        578    1          1   2.8314E-06   2.8314E-06   1.0000E+00  ss.Eqn(349): 0 = v452-(i223)
       503        579    1          1  -1.1710E-09  -1.1710E-09   1.0000E+00  ss.Eqn(350): 0 = v453-(i224)
       504        580    1          1  -1.2694E-06  -1.2694E-06   1.0000E+00  ss.Eqn(351): 0 = v454-(i225)
       505        581    1          1   3.3784E-06   3.3784E-06   1.0000E+00  ss.Eqn(352): 0 = v455-(i226)
       506        582    1          1  -9.4846E-07  -9.4846E-07   1.0000E+00  ss.Eqn(353): 0 = v456-(i227)

A more detailed report of each equation is in the file rto_4_eqn_var.txt that includes a report on each variable that is in the equation. Inequality constraints are converted to equality constraints with an additional slack variable so the equations in gk0_model.apm may look different from those that are written in the Python file. Some of these equations are also from objects such as the m.if3() function.