Issue with sensitivity report using Gekko

Question

I used Gekko for solving a simple linear programming problem. I was able to get an optimal solution. When I tried to get the sensitivity report, I am facing the following error.

@error: Degrees of Freedom Non-Zero for Sensitivity Analysis

However I dont see any problem when I try the same in excel. Please help me find why exactly the issue comes up.

Problem statement :

from gekko import GEKKO
model = GEKKO(remote=False)
model.options.SENSITIVITY = 1   # sensitivity analysis
model.options.SOLVER = 1        # change solver (1=APOPT, 3=IPOPT)
#Maximum demand and implicit contraints as upper bound and lower bound
x1 = model.Var(lb=0, ub=50) # Product DH
x2 = model.Var(lb=0, ub=20) # Product TH
#Objective function
model.Maximize(45*x1+50*x2) # Profit function
#Constraints
model.Equation(500*x1+500*x2<=20000) #Cornflour
model.Equation(750*x1+625*x2<=42000) #Sugar
model.Equation(150*x1+100*x2<=10400) #Fruit and Nut
model.Equation(200*x1+300*x2<=9600) #Ghee
model.solve(disp=True, debug = 1)
p1 = x1.value[0]; p2 = x2.value[0]
print ('Product 1 (DH) in Kgs: ' + str(p1))
print ('Product 2 (TH) in Kgs: ' + str(p2))
print ('Profit       : Rs.' + str(45*p1+50*p2))
print(model.path)

Here is the code output with the error:

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.0489 sec
 Objective      :  -1880.
 Successful solution
 ---------------------------------------------------
 
 
 Generating Sensitivity Analysis
 
 Number of state variables:    6
 Number of total equations: -  4
 Number of independent FVs: -  0
 Number of independent MVs: -  0
 ---------------------------------------
 Degrees of freedom       :    2
 
 Error: DOF must be zero for sensitivity analysis
 @error: Degrees of Freedom Non-Zero for Sensitivity Analysis
 Writing file sensitivity_dof.t0
 STOPPING...

Answer 1

Gekko requires the same number of variables and constraints to perform a sensitivity analysis. It is only for square systems where you are trying to understand the relationship between the inputs to the model and the outputs. Perhaps you need the shadow price of the constraints instead of the sensitivity analysis. Because your problem has only two variables, you can visualize the constraints and objective with a contour plot.

from gekko import GEKKO
model = GEKKO(remote=False)
model.options.SENSITIVITY = 0   # sensitivity analysis
model.options.SOLVER = 1        # change solver (1=APOPT,3=IPOPT)
#Maximum demand and implicit contraints as upper bound and lower bound
x1 = model.Var(lb=0, ub=50) # Product DH
x2 = model.Var(lb=0, ub=20) # Product TH
#Objective function
model.Maximize(45*x1+50*x2) # Profit function
#Constraints
model.Equation(500*x1+500*x2<=20000) #Cornflour
model.Equation(750*x1+625*x2<=42000) #Sugar
model.Equation(150*x1+100*x2<=10400) #Fruit and Nut
model.Equation(200*x1+300*x2<=9600) #Ghee
model.solve(disp=True, debug = 1)
p1 = x1.value[0]; p2 = x2.value[0]
print ('Product 1 (DH) in Kgs: ' + str(p1))
print ('Product 2 (TH) in Kgs: ' + str(p2))
print ('Profit       : Rs.' + str(45*p1+50*p2))
print(model.path)


## Generate a contour plot
# Import some other libraries that we'll need
# matplotlib and numpy packages must also be installed
import matplotlib
import numpy as np
import matplotlib.pyplot as plt

# Design variables at mesh points
x1 = np.arange(0.0, 81.0, 2.0)
x2 = np.arange(0.0, 41.0, 2.0)
x1, x2 = np.meshgrid(x1,x2)

# Equations and Constraints
cf = 500*x1+500*x2
sg = 750*x1+625*x2
fn = 150*x1+100*x2
gh = 200*x1+300*x2

# Objective
obj = 45*x1+50*x2

# Create a contour plot
plt.figure()
# Weight contours
CS = plt.contour(x1,x2,obj)
plt.clabel(CS, inline=1, fontsize=10)
# Constraints
CS = plt.contour(x1, x2, cf,[20000],colors='k',linewidths=[4.0])
plt.clabel(CS, inline=1, fontsize=10)
CS = plt.contour(x1, x2, sg,[42000],colors='b',linewidths=[4.0])
plt.clabel(CS, inline=1, fontsize=10)
CS = plt.contour(x1, x2, fn,[10400],colors='r',linewidths=[4.0])
plt.clabel(CS, inline=1, fontsize=10)
CS = plt.contour(x1, x2, gh,[9600],colors='g',linewidths=[4.0])
plt.clabel(CS, inline=1, fontsize=10)

# plot optimal solution
plt.plot(p1,p2,'o',color='orange',markersize=20)

# plot bounds
plt.plot([0,80],[20,20],'k:',lw=3)
plt.plot([50,50],[0,40],'k:',lw=3)

# Add some labels
plt.title('Linear Programming')
plt.xlabel('Product 1 DH (kgs)')
plt.ylabel('Product 2 TH (kgs)')
plt.grid()
plt.savefig('contour.png')
plt.show()

Edit: Retrieve Shadow Prices (Lagrange Multipliers)

Lagrange multipliers (shadow prices) are available for the constraints when you use the IPOPT solver. Set m.options.SOLVER=3 for IPOPT and the diagnostic level to 2+ with m.options.DIAGLEVEL=2. Here is the modified code to produce the shadow prices:

from gekko import GEKKO
model = GEKKO(remote=False)
model.options.DIAGLEVEL = 2
model.options.SOLVER = 3        # change solver (1=APOPT,3=IPOPT)
#Maximum demand and implicit contraints as upper bound and lower bound
x1 = model.Var(lb=0, ub=50) # Product DH
x2 = model.Var(lb=0, ub=20) # Product TH
#Objective function
model.Maximize(45*x1+50*x2) # Profit function
#Constraints
model.Equation(500*x1+500*x2<=20000) #Cornflour
model.Equation(750*x1+625*x2<=42000) #Sugar
model.Equation(150*x1+100*x2<=10400) #Fruit and Nut
model.Equation(200*x1+300*x2<=9600) #Ghee
model.solve(disp=True, debug = 1)
p1 = x1.value[0]; p2 = x2.value[0]
print ('Product 1 (DH) in Kgs: ' + str(p1))
print ('Product 2 (TH) in Kgs: ' + str(p2))
print ('Profit       : Rs.' + str(45*p1+50*p2))

# for shadow prices, turn on DIAGLEVEL to 2+
#   and use IPOPT solver (APOPT doesn't export Lagrange multipliers)

# Option 1: open the run folder and open apm_lam.txt
model.open_folder()

# Option 2: read apm_lam.txt into array
import numpy as np
lam = np.loadtxt(model.path+'/apm_lam.txt')
print(lam)

The result is:

Product 1 (DH) in Kgs: 24.0
Product 2 (TH) in Kgs: 16.0
Profit       : Rs.1880.0
[-6.99999e-02 -3.44580e-11 -1.01262e-10 -5.00000e-02]