使用 Gekko 的敏感性报告问题
Issue with sensitivity report using Gekko
我使用 Gekko 解决了一个简单的线性规划问题。我能够获得最佳解决方案。当我尝试获取敏感度报告时，我遇到了以下错误。
@error: Degrees of Freedom Non-Zero for Sensitivity Analysis
然而，当我在 excel 中尝试同样的操作时，我没有发现任何问题。
请帮我找出问题出现的确切原因。
问题陈述：
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)
这是错误的代码输出：
---------------------------------------------------
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...
Gekko 需要相同数量的变量和约束来执行敏感性分析。它仅适用于您试图了解模型输入与输出之间关系的方形系统。也许您需要约束的影子价格而不是敏感性分析。因为您的问题只有两个变量，所以您可以使用等高线图可视化约束和 objective。
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()
编辑：检索影子价格（拉格朗日乘数）
当您使用 IPOPT 求解器时，拉格朗日乘数（影子价格）可用于约束。为 IPOPT 设置 m.options.SOLVER=3，并使用 m.options.DIAGLEVEL=2 将诊断级别设置为 2+。这是生成影子价格的修改代码：
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)
结果是：
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]
我使用 Gekko 解决了一个简单的线性规划问题。我能够获得最佳解决方案。当我尝试获取敏感度报告时，我遇到了以下错误。
@error: Degrees of Freedom Non-Zero for Sensitivity Analysis
然而，当我在 excel 中尝试同样的操作时，我没有发现任何问题。 请帮我找出问题出现的确切原因。
问题陈述：
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)
这是错误的代码输出：
---------------------------------------------------
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...
Gekko 需要相同数量的变量和约束来执行敏感性分析。它仅适用于您试图了解模型输入与输出之间关系的方形系统。也许您需要约束的影子价格而不是敏感性分析。因为您的问题只有两个变量，所以您可以使用等高线图可视化约束和 objective。
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()
编辑：检索影子价格（拉格朗日乘数）
当您使用 IPOPT 求解器时，拉格朗日乘数（影子价格）可用于约束。为 IPOPT 设置 m.options.SOLVER=3，并使用 m.options.DIAGLEVEL=2 将诊断级别设置为 2+。这是生成影子价格的修改代码：
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)
结果是：
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]