BYU-PRISM / GEKKO

GEKKO Python for Machine Learning and Dynamic Optimization
https://machinelearning.byu.edu
Other
569 stars 102 forks source link

Add reduced cost, slack or surplus and dual price options to analyze the problem #156

Closed leo-smi closed 1 year ago

leo-smi commented 1 year ago

Like in some no open source tools (LINGO) we have some another options to analze our model results:

image

APMonitor commented 1 year ago

Here is information on how to report Lagrange multipliers (Dual Price): https://stackoverflow.com/questions/62040814/lagrange-multipliers-marginals-in-gekko The solvers in Gekko are all nonlinear so the Dual price may also be nonlinear. The Lagrange multiplier is a linear approximation of shadow (dual) price.

Slack (surplus) variables are created with each inequality constraint. They are found in the solution summary report in m.path in the .t0 file for that mode (IMODE=3 is rto.t0).

I appreciate the suggestion to create a summary. I've added it as a feature request.

leo-smi commented 1 year ago

Thank you so much

leo-smi commented 1 year ago

Sorry to bother you, could you do an example to obtain the .t0 file? I got the apm_lam.txt with the Lagrange values, but I can't load the .t0 (I have set the prob.options.IMODE = 3)

Here's my code:

from gekko import GEKKO
import numpy as np

# Define environment

prob = GEKKO(remote=False)
prob.options.IMODE = 3
print(prob.path)

# Define decision variables

A = prob.Var(lb=0, ub=None, integer=True) 
B = prob.Var(lb=0, ub=None, integer=True)
C = prob.Var(lb=0, ub=None, integer=True)

# Add objective function to the environment
# m.Obj(obj)
# m.Minimize(-obj) == m.Maximize(obj)

obj_func = (200*A + 300*B + 400*C)
prob.Obj(obj_func) # Objective by default is same as prob.Minimize(Obj)
prob.options.DIAGLEVEL=2 # The Lagrange multiplier is a linear approximation of shadow (dual) price

# Add constraints to the environment

constraints = [
    CAM_A :=       A                       <= 4,
    CAM_B :=                  B            <= 4,
    CAM_C :=                             C <= 2,
    CARGA := 5_000*A + 10_000*B + 20_000*C >= 80_000,
]

for const in constraints:
    prob.Equation(const)

print(prob.Equation)   

# Solve the problem (1: MINLP solver, 2,3: Other Solvers)

prob.options.SOLVER=1
prob.solve(disp=False)

# Display results

print('Results')
print(f'A: {A.value}')
print(f'B: {B.value}')
print(f'C: {C.value}')
print(f'Optimal Value of Objective Is = {prob.options.objfcnval}')

# https://github.com/BYU-PRISM/GEKKO/issues/156#issuecomment-1365980918
# The Lagrange multiplier is a linear approximation of shadow (dual) price
print('Lagrange multipliers')
# lam = np.loadtxt(prob.path + '/apm_lam.txt')
lam = np.loadtxt(prob.path + '/rto.t0')
print(lam)

I got some error:

# https://github.com/BYU-PRISM/GEKKO/issues/156#issuecomment-1365980918
# The Lagrange multiplier is a linear approximation of shadow (dual) price
print('Lagrange multipliers')
# lam = np.loadtxt(prob.path + '/apm_lam.txt')
lam = np.loadtxt(prob.path + '/rto.t0')
print(lam)
ValueError                                Traceback (most recent call last)
Cell In [193], line 5
      3 print('Lagrange multipliers')
      4 # lam = np.loadtxt(prob.path + '/apm_lam.txt')
----> 5 lam = np.loadtxt(prob.path + '/rto.t0')
      6 print(lam)

ValueError: could not convert string 'Value' to float64 at row 0, column 1.

Acctually I have no idead what I'm doing.

APMonitor commented 1 year ago

Here is a version that reports the slack variables:

from gekko import GEKKO
import numpy as np
import pandas as pd

# Define environment

prob = GEKKO(remote=False)
prob.options.IMODE = 3
print(prob.path)

# Define decision variables

A = prob.Var(lb=0, ub=None, integer=True, name='A') 
B = prob.Var(lb=0, ub=None, integer=True, name='B')
C = prob.Var(lb=0, ub=None, integer=True, name='C')

# Add objective function to the environment
# m.Obj(obj)
# m.Minimize(-obj) == m.Maximize(obj)

obj_func = (200*A + 300*B + 400*C)
prob.Minimize(obj_func) # Objective by default is same as prob.Minimize(Obj)
prob.options.DIAGLEVEL=2 # The Lagrange multiplier is a linear approximation of shadow (dual) price

# Add constraints to the environment

constraints = [
    CAM_A :=       A                       <= 4,
    CAM_B :=                  B            <= 4,
    CAM_C :=                             C <= 2,
    CARGA := 5_000*A + 10_000*B + 20_000*C >= 80_000,
]

for const in constraints:
    prob.Equation(const)

print(prob.Equation)   

# Solve the problem (1: MINLP solver, 2,3: Other Solvers)

prob.options.SOLVER=1
prob.solve(disp=False)

# Display results

print('Results')
print(f'A: {A.value[0]}')
print(f'B: {B.value[0]}')
print(f'C: {C.value[0]}')
print(f'Optimal Value of Objective Is = {prob.options.objfcnval}')

# https://github.com/BYU-PRISM/GEKKO/issues/156#issuecomment-1365980918
# The Lagrange multiplier is a linear approximation of shadow (dual) price
print('Lagrange multipliers')
lam = np.loadtxt(prob.path + '/apm_lam.txt')
print(lam)
results = pd.read_csv(prob.path + '/results.csv',header=None)
results.columns=['Name','Value']
print(results)

The result is shown below. The results.csv file is easier to parse, although rto.t0 has the same information. The int_ in front of the variables A, B, and C indicates to APMonitor that it should be treated as an Integer type decision variable. The slk_ in front of the last four variables indicates that they are slack variables with lower bound of 0. These correspond to the four inequality equations defined in the code.

Results
A: 0.0
B: 4.0
C: 2.0
Optimal Value of Objective Is = 2000.0
Lagrange multipliers
[0. 0. 0. 0.]
    Name  Value
0  int_a    0.0
1  int_b    4.0
2  int_c    2.0
3  slk_1    4.0
4  slk_2    0.0
5  slk_3    0.0
6  slk_4    0.0

Some of the solvers, such as APOPT, don't report Lagrange multipliers. Solving with IPOPT (prob.options.SOLVER=3) gives the correct value of the Lagrange multipliers. It is nearly an integer solution although IPOPT is an NLP solver.

Results
A: 0.0
B: 3.9999999842
C: 2.0000000096
Optimal Value of Objective Is = 1999.9999991
Lagrange multipliers
[-1.23654978e-06 -3.67116766e-01 -5.07342333e+01 -7.53671155e-03]
    Name         Value
0  int_a  0.000000e+00
1  int_b  4.000000e+00
2  int_c  2.000000e+00
3  slk_1  4.000000e+00
4  slk_2  1.579791e-08
5  slk_3  0.000000e+00
6  slk_4  2.973699e-06
leo-smi commented 1 year ago

Thank you so much, the results match.

image