标题可能措辞不当。我需要的是根据分析中包含的“国家”数量来扩展或收缩的编码。该问题至少需要 2 个国家，但根据所解决的问题可以包括 N 个国家。例如，使用 scipy.optimize，我可以为 2 个国家或 5 个国家/地区进行编码，但我无法弄清楚如何编写代码，以便同一组指令可以在第一遍中容纳 2 个国家，然后在第二遍中容纳 5 个国家/地区。如果我的所有方程都是线性的，我可以轻松地在 X 矩阵和 Y 向量中添加或减去，以求解任意数量国家的 np.linalg.inv(X).dot(Y) 。我想使用 scipy.optimize 做同样的事情，因为我的方程是非线性的。请注意，我对 scipy.optimize 并不感兴趣。我会使用任何模块，只要它能解决问题，并且只需提供 N 个国家的数据即可扩展或缩小问题。
以下编码工作得很好。但是假设我想要 Trade=K=6 而不是 4。我该如何编写这个程序，以便可以在不重写程序的情况下求解 K = 4 或 6？
还有，我是一个通过Fortran课程学习编程的老家伙。我认为我的 python 编码看起来很笨重。我会接受任何关于如何使这段代码更优雅的建议。
import numpy as np
from scipy.optimize import minimize
Trade = int(4)
K = Trade
P0 = np.zeros(K+1)
Q0 = np.zeros(K+1)
Qc0 = np.zeros(K+1)
Qex0 = np.zeros(K)
TSup = np.zeros(K)
Qp0 = np.zeros(K+1)
Qim0 = np.zeros(K)
ED = np.zeros(K+1)
ES = np.zeros(K+1)
P0[0] = float( 81.8729) # Initial Price U.S.
P0[1] = float( 51.6940) # Initial Price Mexico
P0[2] = float( 96.9458) # Initial Price Canada
P0[3] = float( 61.7404) # Initial Price Rest of the World (ROW)
# Note that Quantity Consumed is equal to Quantity Demanded
Qc0[0] = float( 99861.9890) # Quantity Consumed U.S.
Qc0[1] = float( 13164.0000) # Quantity Consumed Mexico
Qc0[2] = float( 10490.0000) # Quantity Consumed Canada
Qc0[3] = float( 84311.5110) # Quantity Consumed ROW
Qex0[0] = float( 6060.9012) # Quantity Exported from the U.S.
Qex0[1] = float( 13.0000) # Quantity Exported from Mexico
Qex0[2] = float( 20.0000) # Quantity Exported from Canada
Qex0[3] = float( 2537515.7206) # Quantity Exported from the ROW
# Note that Total Supply is equal to Total Distributions
TSup[0] = float( 105922.8902) # Total Supply U.S.
TSup[1] = float( 13177.0000) # Total Supply Mexico
TSup[2] = float( 10530.0000) # Total Supply Canada
TSup[3] = float( 2621807.2316) # Total Supply ROW
# Note that Quantity Produces is equal to Quantity Supplied
Qp0[0] = float( 102721.4352) # Quantity Produced U.S.
Qp0[1] = float( 13152.0000) # Quantity Produced Mexico
Qp0[2] = float( 10440.0000) # Quantity Produced Canada
Qp0[3] = float( 420501.0646) # Quantity Produced ROW
Qim0[0] = float( 3201.4550) # Quantity Imported into the U.S.
Qim0[1] = float( 25.0000) # Quantity Imported into Mexico
Qim0[2] = float( 90.0000) # Quantity Imported into Canada
Qim0[3] = float( 2201306.1670) # Quantity Imported into the ROW
ED[0] = float( -1.65 ) # U.S. ED for Whole Milk E.U. All Dairy -0.57
ED[1] = float( -0.68 ) # Mexico ED for Fluid Milk
ED[2] = float( -0.5125) # Canada ED for Milk Products
ED[3] = float( -1.1283) # ROW ED
ES[0] = float( 2.25 ) # U.S. ES for Milk
ES[1] = float( 1.128343478)
ES[2] = float( 1.689 )
ES[3] = float( 1.689 )
DDP = np.zeros(K) # A Direct Change in Demand Price
DDQc = np.zeros(K) # A Direct Change in Quantity Consumed (e.g. Income Effects)
DDQim = np.zeros(K) # A Direct Change in Quantity Imported
DSP = np.zeros(K) # A Direct Change in Supply Price
DSQp = np.zeros(K) # A Direct Change in Quantity Produced
#DSQ[0] = float( 27750.3412)
DSQex = np.zeros(K) # A Direct Change in Quantity Exported
### Start Non-linear Optimization
### to Determine Initial Values
###################################
###################################
###################################
Xx = np.zeros(6+5*K)
# The Objective Function: Maximize World Surplus
def objective(x):
return -1*((((Xx[1] -Xx[2]) *Xx[5]) /2)+ \
(((Xx[6] -Xx[7]) *Xx[10])/2)+ \
(((Xx[11]-Xx[12])*Xx[15])/2)+ \
(((Xx[16]-Xx[17])*Xx[20])/2))
### Start Calculation of Starting Values (initial guesses)
###################################
###################################
P0[K] = np.sum(np.delete(P0, [K]))/K
Qc0[K] = float( 1.0000)
Qp0[K] = float( 1.0000)
Q0 = (Qc0+Qp0)/2
Q0[K] = np.sum(np.delete(Q0, [K]))
# Starting Values for Demand
ds = np.zeros(K+1)
ED[K] = float(-1.0000)
ds = (1/ED)*(P0/Q0)
dc = np.zeros(K)
dc = P0 - ds*Q0
dc[K] = np.max(np.delete(dc, [K]))
ds[K] = (P0[K]-dc[K])/Q0[K]
# Starting Values for Supply
ss = np.zeros(K+1)
ES[K] = float( 1.0000)
ss = (1/ES)*(P0/Q0)
sc = np.zeros(K)
sc = P0 - ss*Q0
sc[K] = np.min(np.delete(sc, [K]))
ss[K] = (P0[K]-sc[K])/Q0[K]
SV = np.zeros(6+5*K)
k = int(0)
for k in range(K+1):
if k == 0: SV[k]=P0[K]
SV[1+5*k] = dc[k]
SV[2+5*k] = sc[k]
SV[3+5*k] = ds[k]
SV[4+5*k] = ss[k]
SV[5+5*k] = Q0[k]
###################################
###################################
### End Calculation of Starting Values
# Define variable lower and upper bounds
bnds = ((0.0,None),
(0.0,None),
(None,None),
(None,0.0),
(0.0,None),
(0.0,None),
(0.0,None),
(None,None),
(None,0.0),
(0.0,None),
(0.0,None),
(0.0,None),
(None,None),
(None,0.0),
(0.0,None),
(0.0,None),
(0.0,None),
(None,None),
(None,0.0),
(0.0,None),
(0.0,None),
(0.0,None),
(None,None),
(None,0.0),
(0.0,None),
(0.0,None))
### Start Constraint Work
###################################
###################################
# Define the Constraints for Demand
#Con = ["Con"+str(x+1)+"(x)" for x in range(K)]
def Con000(x):
return Xx[1] +Xx[3] *Xx[5] -P0[0]
def Con001(x):
return Xx[0]-Xx[1] -Xx[3] *Qc0[0]
def Con002(x):
return Xx[6] +Xx[8] *Xx[10]-P0[1]
def Con003(x):
return Xx[0]-Xx[6] -Xx[8] *Qc0[1]
def Con004(x):
return Xx[11]+Xx[13]*Xx[15]-P0[2]
def Con005(x):
return Xx[0]-Xx[11]-Xx[13]*Qc0[2]
def Con006(x):
return Xx[16]+Xx[18]*Xx[20]-P0[3]
def Con007(x):
return Xx[0]-Xx[16]-Xx[18]*Qc0[3]
def Con008(x):
return Xx[0]-Xx[21]-Xx[23]*Xx[25]
def Con010(x):
return Xx[21]-(Xx[1]+Xx[6]+Xx[11]+Xx[16])
def Con011(x):
return Xx[25]- (Xx[5]+Xx[10]+Xx[15]+Xx[20])
def Con020(x):
return P0[0]-Xx[1]
def Con021(x):
return P0[1]-Xx[6]
def Con022(x):
return P0[2]-Xx[11]
def Con023(x):
return P0[3]-Xx[16]
def Con024(x):
return Xx[0]-Xx[1]
def Con025(x):
return Xx[0]-Xx[6]
def Con026(x):
return Xx[0]-Xx[11]
def Con027(x):
return Xx[0]-Xx[16]
def Con028(x):
return Xx[0]-Xx[21]
# Define the Constraints for Supply
def Con100(x):
return Xx[2] +Xx[4] *Xx[5] -P0[0]
def Con101(x):
return Xx[0]-Xx[2] -Xx[4] *Qp0[0]
def Con102(x):
return Xx[7] +Xx[9] *Xx[10]-P0[1]
def Con103(x):
return Xx[0]-Xx[7] -Xx[9] *Qp0[1]
def Con104(x):
return Xx[12]+Xx[14]*Xx[15]-P0[2]
def Con105(x):
return Xx[0]-Xx[12]-Xx[14]*Qp0[2]
def Con106(x):
return Xx[17]+Xx[19]*Xx[20]-P0[3]
def Con107(x):
return Xx[0]-Xx[17]-Xx[19]*Qp0[3]
def Con108(x):
return Xx[0]-Xx[22]-Xx[24]*Xx[25]
def Con110(x):
return Xx[22]-(Xx[2]+Xx[7]+Xx[12]+Xx[17])
def Con120(x):
return Xx[2] -P0[0]
def Con121(x):
return Xx[7] -P0[1]
def Con122(x):
return Xx[12]-P0[2]
def Con123(x):
return Xx[17]-P0[3]
def Con124(x):
return Xx[2] -Xx[0]
def Con125(x):
return Xx[7] -Xx[0]
def Con126(x):
return Xx[12]-Xx[0]
def Con127(x):
return Xx[17]-Xx[0]
def Con128(x):
return Xx[22]-Xx[0]
if Qp0[0]-Qc0[0] >= 0:
def Con030(x):
return Qc0[0]-Xx[5]
def Con130(x):
return Xx[5]-Qp0[0]
else:
def Con030(x):
return Xx[5]-Qc0[0]
def Con130(x):
return Qp0[0]-Xx[5]
if Qp0[1]-Qc0[1] >= 0:
def Con031(x):
return Qc0[1]-Xx[10]
def Con131(x):
return Xx[10]-Qp0[1]
else:
def Con031(x):
return Xx[10]-Qc0[1]
def Con131(x):
return Qp0[1]-Xx[10]
if Qp0[2]-Qc0[2] >= 0:
def Con032(x):
return Qc0[2]-Xx[15]
def Con132(x):
return Xx[15]-Qp0[2]
else:
def Con032(x):
return Xx[15]-Qc0[2]
def Con132(x):
return Qp0[2]-Xx[15]
if Qp0[3]-Qc0[3] >= 0:
def Con033(x):
return Qc0[3]-Xx[20]
def Con133(x):
return Xx[20]-Qp0[3]
else:
def Con033(x):
return Xx[20]-Qc0[3]
def Con133(x):
return Qp0[3]-Xx[20]
if np.sum(Qp0)-np.sum(Qc0) >= 0:
def Con034(x):
return np.sum(Qc0)-Xx[25]
def Con134(x):
return Xx[25]-np.sum(Qp0)
else:
def Con034(x):
return Xx[25]-np.sum(Qc0)
def Con134(x):
return np.sum(Qp0)-Xx[25]
def Con200(x):
return -1*((((Xx[1] -Xx[2]) *Xx[5]) /2) + \
(((Xx[6] -Xx[7]) *Xx[10])/2) + \
(((Xx[11]-Xx[12])*Xx[15])/2) + \
(((Xx[16]-Xx[17])*Xx[20])/2))+ \
(((Xx[21]-Xx[22])*Xx[25])/2)
# Type and Collate the Constraints
C_D11 = {'type': 'eq', 'fun': Con000}
C_D12 = {'type': 'eq', 'fun': Con001}
C_D13 = {'type': 'eq', 'fun': Con002}
C_D14 = {'type': 'eq', 'fun': Con003}
C_D15 = {'type': 'eq', 'fun': Con004}
C_D16 = {'type': 'eq', 'fun': Con005}
C_D17 = {'type': 'eq', 'fun': Con006}
C_D18 = {'type': 'eq', 'fun': Con007}
C_D19 = {'type': 'eq', 'fun': Con008}
C_D21 = {'type': 'eq', 'fun': Con010}
C_D22 = {'type': 'eq', 'fun': Con011}
C_D31 = {'type': 'ineq', 'fun': Con020}
C_D32 = {'type': 'ineq', 'fun': Con021}
C_D33 = {'type': 'ineq', 'fun': Con022}
C_D34 = {'type': 'ineq', 'fun': Con023}
C_D35 = {'type': 'ineq', 'fun': Con024}
C_D36 = {'type': 'ineq', 'fun': Con025}
C_D37 = {'type': 'ineq', 'fun': Con026}
C_D38 = {'type': 'ineq', 'fun': Con027}
C_D39 = {'type': 'ineq', 'fun': Con028}
C_D41 = {'type': 'ineq', 'fun': Con030}
C_D42 = {'type': 'ineq', 'fun': Con031}
C_D43 = {'type': 'ineq', 'fun': Con032}
C_D44 = {'type': 'ineq', 'fun': Con033}
C_D45 = {'type': 'ineq', 'fun': Con034}
C_S11 = {'type': 'eq', 'fun': Con100}
C_S12 = {'type': 'eq', 'fun': Con101}
C_S13 = {'type': 'eq', 'fun': Con102}
C_S14 = {'type': 'eq', 'fun': Con103}
C_S15 = {'type': 'eq', 'fun': Con104}
C_S16 = {'type': 'eq', 'fun': Con105}
C_S17 = {'type': 'eq', 'fun': Con106}
C_S18 = {'type': 'eq', 'fun': Con107}
C_S19 = {'type': 'eq', 'fun': Con108}
C_S21 = {'type': 'eq', 'fun': Con110}
C_S31 = {'type': 'ineq', 'fun': Con120}
C_S32 = {'type': 'ineq', 'fun': Con121}
C_S33 = {'type': 'ineq', 'fun': Con122}
C_S34 = {'type': 'ineq', 'fun': Con123}
C_S35 = {'type': 'ineq', 'fun': Con124}
C_S36 = {'type': 'ineq', 'fun': Con125}
C_S37 = {'type': 'ineq', 'fun': Con126}
C_S38 = {'type': 'ineq', 'fun': Con127}
C_S39 = {'type': 'ineq', 'fun': Con128}
C_S41 = {'type': 'ineq', 'fun': Con130}
C_S42 = {'type': 'ineq', 'fun': Con131}
C_S43 = {'type': 'ineq', 'fun': Con132}
C_S44 = {'type': 'ineq', 'fun': Con133}
C_S45 = {'type': 'ineq', 'fun': Con134}
C_S46 = {'type': 'eq', 'fun': Con200}
Cons = ([C_D11, C_D12, C_D13, C_D14, C_D15, C_D16, C_D17, C_D18, C_D19, \
C_D21, C_D22, \
C_D31, C_D32, C_D33, C_D34, C_D35, C_D36, C_D37, C_D38, C_D39, \
C_D41, C_D42, C_D43, C_D44, C_D45,
C_S11, C_S12, C_S13, C_S14, C_S15, C_S16, C_S17, C_S18, C_S19, \
C_S21, \
C_S31, C_S32, C_S33, C_S34, C_S35, C_S36, C_S37, C_S38, C_S39, \
C_S41, C_S42, C_S43, C_S44, C_S45])
###################################
###################################
### End Constraint Work
solution = minimize(objective,SV,method='SLSQP',\
bounds=bnds,constraints=Cons)
Xx = solution.x
# Return Model Results to Variables
k = int(0)
for k in range(K):
if k == 0: P0[k] = Xx[k]
dc[k] = Xx[1+5*k]
sc[k] = Xx[2+5*k]
ds[k] = Xx[3+5*k]
ss[k] = Xx[4+5*k]
Q0[k] = Xx[5+5*k]
print ("Initial Prices")
print (P0)
print (' ')
print ("Initial Demand Choke Price")
print (dc)
print ("Initial Supply Choke Prices")
print (sc)
print (' ')
print ("Initial Slope of the Demand Curve")
print (ds)
print (' ')
print ("Initial Slope of the Supply Curve")
print (ss)
print (' ')
print ("Initial Quantity")
print (Q0)
print (' ')
Qc0 = np.delete(Qc0, [K])
print ("Initial Quantity Consumed")
print (Qc0)
print (' ')
Qp0 = np.delete(Qp0, [K])
print ("Initial Quantity Produced")
print (Qp0)
print (' ')
dsT = np.zeros(K+1)
ED[K] = float(-1.0000)
dsT = (1/ED)*(P0/Q0)
ds[K] = (P0[K]-dc[K])/Q0[K]
print (' ')
print ("Initial Slope of the Demand Curve")
print (ds)
print ("Test Slope of the Demand Curve")
print (dsT)
print (' ')
ssT = np.zeros(K+1)
ES[K] = float( 1.0000)
ssT = (1/ES)*(P0/Q0)
ss[K] = (P0[K]-sc[K])/Q0[K]
print (' ')
print ("Initial Slope of the Supply Curve")
print (ss)
print ("Test Slope of the Supply Curve")
print (ssT)
print (' ')
###################################
###################################
###################################
### End Non-linear Optimization
我将停止完整的演示，只写一个需要更改的列表。完成繁重的工作，如果您有一些不起作用的问题，请发布一个新问题（并将其链接到此处）；如果您有一些确实可以工作的内容，请改为在 codereview.stackexchange.com 上发布。
float(81.8)81.8为什么要向一半向量添加一个元素？如果它确实是“起始值”，那么它不属于您的其他数据。它应该使用一次，将
x0minimize摆脱顶部输入数据的所有常量索引分配。编写此 CSV 或类似的 CSV：
country, initial price, consumed, exported, supply, produced, imported, ED, ES
US, 81.8729, 99861.9890, 6060.9012, 105922.8902, 102721.4352, 3201.4550, -1.6500, 2.250000000
Mexico, 51.6940, 13164.0000, 13.0000, 13177.0000, 13152.0000, 25.0000, -0.6800, 1.128343478
Canada, 96.9458, 10490.0000, 20.0000, 10530.0000, 10440.0000, 90.0000, -0.5125, 1.689000000
ROW, 61.7404, 84311.5110, 2537515.7206, 2621807.2316, 420501.0646, 2201306.1670, -1.1283, 1.689000000
然后使用
pd.read_csv你的目标和每一项限制都是错误的。特别注意
objective(x)(x)x，而不是 Xx实际决策变量。为了强制执行这一点并阻止意外的闭包绑定，我强烈建议您将所有代码移至没有全局状态的函数。
用切片
np.delete(P0, [K])P0[:-1]LinearConstraintNonlinearConstraint不要在括号表达式内部写续行
\添加回归测试。
执行上述一小部分操作的编辑可能看起来像这样
import numpy as np
from scipy.optimize import OptimizeResult
def main() -> None:
Trade = int(4)
K = Trade
P0 = np.zeros(K + 1)
Qc0 = np.zeros(K + 1)
Qex0 = np.zeros(K)
TSup = np.zeros(K)
Qp0 = np.zeros(K + 1)
Qim0 = np.zeros(K)
ED = np.zeros(K + 1)
ES = np.zeros(K + 1)
P0[0] = 81.8729 # Initial Price U.S.
P0[1] = 51.6940 # Initial Price Mexico
P0[2] = 96.9458 # Initial Price Canada
P0[3] = 61.7404 # Initial Price Rest of the World (ROW)
# Note that Quantity Consumed is equal to Quantity Demanded
Qc0[0] = 99861.9890 # Quantity Consumed U.S.
Qc0[1] = 13164.0000 # Quantity Consumed Mexico
Qc0[2] = 10490.0000 # Quantity Consumed Canada
Qc0[3] = 84311.5110 # Quantity Consumed ROW
Qex0[0] = 6060.9012 # Quantity Exported from the U.S.
Qex0[1] = 13.0000 # Quantity Exported from Mexico
Qex0[2] = 20.0000 # Quantity Exported from Canada
Qex0[3] = 2537515.7206 # Quantity Exported from the ROW
# Note that Total Supply is equal to Total Distributions
TSup[0] = 105922.8902 # Total Supply U.S.
TSup[1] = 13177.0000 # Total Supply Mexico
TSup[2] = 10530.0000 # Total Supply Canada
TSup[3] = 2621807.2316 # Total Supply ROW
# Note that Quantity Produces is equal to Quantity Supplied
Qp0[0] = 102721.4352 # Quantity Produced U.S.
Qp0[1] = 13152.0000 # Quantity Produced Mexico
Qp0[2] = 10440.0000 # Quantity Produced Canada
Qp0[3] = 420501.0646 # Quantity Produced ROW
Qim0[0] = 3201.4550 # Quantity Imported into the U.S.
Qim0[1] = 25.0000 # Quantity Imported into Mexico
Qim0[2] = 90.0000 # Quantity Imported into Canada
Qim0[3] = 2201306.1670 # Quantity Imported into the ROW
ED[0] = -1.65 # U.S. ED for Whole Milk E.U. All Dairy -0.57
ED[1] = -0.68 # Mexico ED for Fluid Milk
ED[2] = -0.5125 # Canada ED for Milk Products
ED[3] = -1.1283 # ROW ED
ES[0] = 2.25 # U.S. ES for Milk
ES[1] = 1.128343478
ES[2] = 1.689
ES[3] = 1.689
Xx = np.zeros(6 + 5*K)
# The Objective Function: Maximize World Surplus
def objective(x):
return -1 * ((((Xx[1] - Xx[2]) * Xx[5]) / 2) +
(((Xx[6] - Xx[7]) * Xx[10]) / 2) +
(((Xx[11] - Xx[12]) * Xx[15]) / 2) +
(((Xx[16] - Xx[17]) * Xx[20]) / 2))
### Start Calculation of Starting Values (initial guesses)
P0[K] = np.sum(P0) / K
Qc0[K] = 1.0000
Qp0[K] = 1.0000
Q0 = (Qc0 + Qp0) / 2
Q0[K] = np.sum(Q0)
# Starting Values for Demand
ED[K] = -1.0000
ds = (1 / ED) * (P0 / Q0)
dc = P0 - ds * Q0
dc[K] = np.max(dc)
ds[K] = (P0[K] - dc[K]) / Q0[K]
# Starting Values for Supply
ES[K] = 1.0000
ss = (1 / ES) * (P0 / Q0)
sc = P0 - ss * Q0
sc[K] = np.min(np.delete(sc, [K]))
ss[K] = (P0[K] - sc[K]) / Q0[K]
SV = np.zeros(6 + 5 * K)
for k in range(K + 1):
if k == 0: SV[k] = P0[K]
SV[1 + 5 * k] = dc[k]
SV[2 + 5 * k] = sc[k]
SV[3 + 5 * k] = ds[k]
SV[4 + 5 * k] = ss[k]
SV[5 + 5 * k] = Q0[k]
# Define variable lower and upper bounds
bnds = ((0.0, None),
(0.0, None),
(None, None),
(None, 0.0),
(0.0, None),
(0.0, None),
(0.0, None),
(None, None),
(None, 0.0),
(0.0, None),
(0.0, None),
(0.0, None),
(None, None),
(None, 0.0),
(0.0, None),
(0.0, None),
(0.0, None),
(None, None),
(None, 0.0),
(0.0, None),
(0.0, None),
(0.0, None),
(None, None),
(None, 0.0),
(0.0, None),
(0.0, None))
### Start Constraint Work
# Define the Constraints for Demand
# Con = ["Con"+str(x+1)+"(x)" for x in range(K)]
def Con000(x):
return Xx[1] + Xx[3] * Xx[5] - P0[0]
def Con001(x):
return Xx[0] - Xx[1] - Xx[3] * Qc0[0]
def Con002(x):
return Xx[6] + Xx[8] * Xx[10] - P0[1]
def Con003(x):
return Xx[0] - Xx[6] - Xx[8] * Qc0[1]
def Con004(x):
return Xx[11] + Xx[13] * Xx[15] - P0[2]
def Con005(x):
return Xx[0] - Xx[11] - Xx[13] * Qc0[2]
def Con006(x):
return Xx[16] + Xx[18] * Xx[20] - P0[3]
def Con007(x):
return Xx[0] - Xx[16] - Xx[18] * Qc0[3]
def Con008(x):
return Xx[0] - Xx[21] - Xx[23] * Xx[25]
def Con010(x):
return Xx[21] - (Xx[1] + Xx[6] + Xx[11] + Xx[16])
def Con011(x):
return Xx[25] - (Xx[5] + Xx[10] + Xx[15] + Xx[20])
def Con020(x):
return P0[0] - Xx[1]
def Con021(x):
return P0[1] - Xx[6]
def Con022(x):
return P0[2] - Xx[11]
def Con023(x):
return P0[3] - Xx[16]
def Con024(x):
return Xx[0] - Xx[1]
def Con025(x):
return Xx[0] - Xx[6]
def Con026(x):
return Xx[0] - Xx[11]
def Con027(x):
return Xx[0] - Xx[16]
def Con028(x):
return Xx[0] - Xx[21]
# Define the Constraints for Supply
def Con100(x):
return Xx[2] + Xx[4] * Xx[5] - P0[0]
def Con101(x):
return Xx[0] - Xx[2] - Xx[4] * Qp0[0]
def Con102(x):
return Xx[7] + Xx[9] * Xx[10] - P0[1]
def Con103(x):
return Xx[0] - Xx[7] - Xx[9] * Qp0[1]
def Con104(x):
return Xx[12] + Xx[14] * Xx[15] - P0[2]
def Con105(x):
return Xx[0] - Xx[12] - Xx[14] * Qp0[2]
def Con106(x):
return Xx[17] + Xx[19] * Xx[20] - P0[3]
def Con107(x):
return Xx[0] - Xx[17] - Xx[19] * Qp0[3]
def Con108(x):
return Xx[0] - Xx[22] - Xx[24] * Xx[25]
def Con110(x):
return Xx[22] - (Xx[2] + Xx[7] + Xx[12] + Xx[17])
def Con120(x):
return Xx[2] - P0[0]
def Con121(x):
return Xx[7] - P0[1]
def Con122(x):
return Xx[12] - P0[2]
def Con123(x):
return Xx[17] - P0[3]
def Con124(x):
return Xx[2] - Xx[0]
def Con125(x):
return Xx[7] - Xx[0]
def Con126(x):
return Xx[12] - Xx[0]
def Con127(x):
return Xx[17] - Xx[0]
def Con128(x):
return Xx[22] - Xx[0]
if Qp0[0] - Qc0[0] >= 0:
def Con030(x):
return Qc0[0] - Xx[5]
def Con130(x):
return Xx[5] - Qp0[0]
else:
def Con030(x):
return Xx[5] - Qc0[0]
def Con130(x):
return Qp0[0] - Xx[5]
if Qp0[1] - Qc0[1] >= 0:
def Con031(x):
return Qc0[1] - Xx[10]
def Con131(x):
return Xx[10] - Qp0[1]
else:
def Con031(x):
return Xx[10] - Qc0[1]
def Con131(x):
return Qp0[1] - Xx[10]
if Qp0[2] - Qc0[2] >= 0:
def Con032(x):
return Qc0[2] - Xx[15]
def Con132(x):
return Xx[15] - Qp0[2]
else:
def Con032(x):
return Xx[15] - Qc0[2]
def Con132(x):
return Qp0[2] - Xx[15]
if Qp0[3] - Qc0[3] >= 0:
def Con033(x):
return Qc0[3] - Xx[20]
def Con133(x):
return Xx[20] - Qp0[3]
else:
def Con033(x):
return Xx[20] - Qc0[3]
def Con133(x):
return Qp0[3] - Xx[20]
if np.sum(Qp0) - np.sum(Qc0) >= 0:
def Con034(x):
return np.sum(Qc0) - Xx[25]
def Con134(x):
return Xx[25] - np.sum(Qp0)
else:
def Con034(x):
return Xx[25] - np.sum(Qc0)
def Con134(x):
return np.sum(Qp0) - Xx[25]
def Con200(x):
return (
-(
(Xx[1] - Xx[2]) * Xx[5]/2 +
(Xx[6] - Xx[7]) * Xx[10]/2 +
(Xx[11] - Xx[12]) * Xx[15]/2 +
(Xx[16] - Xx[17]) * Xx[20]/2
) +
(Xx[21] - Xx[22]) * Xx[25] / 2
)
# Type and Collate the Constraints
C_D11 = {'type': 'eq', 'fun': Con000}
C_D12 = {'type': 'eq', 'fun': Con001}
C_D13 = {'type': 'eq', 'fun': Con002}
C_D14 = {'type': 'eq', 'fun': Con003}
C_D15 = {'type': 'eq', 'fun': Con004}
C_D16 = {'type': 'eq', 'fun': Con005}
C_D17 = {'type': 'eq', 'fun': Con006}
C_D18 = {'type': 'eq', 'fun': Con007}
C_D19 = {'type': 'eq', 'fun': Con008}
C_D21 = {'type': 'eq', 'fun': Con010}
C_D22 = {'type': 'eq', 'fun': Con011}
C_D31 = {'type': 'ineq', 'fun': Con020}
C_D32 = {'type': 'ineq', 'fun': Con021}
C_D33 = {'type': 'ineq', 'fun': Con022}
C_D34 = {'type': 'ineq', 'fun': Con023}
C_D35 = {'type': 'ineq', 'fun': Con024}
C_D36 = {'type': 'ineq', 'fun': Con025}
C_D37 = {'type': 'ineq', 'fun': Con026}
C_D38 = {'type': 'ineq', 'fun': Con027}
C_D39 = {'type': 'ineq', 'fun': Con028}
C_D41 = {'type': 'ineq', 'fun': Con030}
C_D42 = {'type': 'ineq', 'fun': Con031}
C_D43 = {'type': 'ineq', 'fun': Con032}
C_D44 = {'type': 'ineq', 'fun': Con033}
C_D45 = {'type': 'ineq', 'fun': Con034}
C_S11 = {'type': 'eq', 'fun': Con100}
C_S12 = {'type': 'eq', 'fun': Con101}
C_S13 = {'type': 'eq', 'fun': Con102}
C_S14 = {'type': 'eq', 'fun': Con103}
C_S15 = {'type': 'eq', 'fun': Con104}
C_S16 = {'type': 'eq', 'fun': Con105}
C_S17 = {'type': 'eq', 'fun': Con106}
C_S18 = {'type': 'eq', 'fun': Con107}
C_S19 = {'type': 'eq', 'fun': Con108}
C_S21 = {'type': 'eq', 'fun': Con110}
C_S31 = {'type': 'ineq', 'fun': Con120}
C_S32 = {'type': 'ineq', 'fun': Con121}
C_S33 = {'type': 'ineq', 'fun': Con122}
C_S34 = {'type': 'ineq', 'fun': Con123}
C_S35 = {'type': 'ineq', 'fun': Con124}
C_S36 = {'type': 'ineq', 'fun': Con125}
C_S37 = {'type': 'ineq', 'fun': Con126}
C_S38 = {'type': 'ineq', 'fun': Con127}
C_S39 = {'type': 'ineq', 'fun': Con128}
C_S41 = {'type': 'ineq', 'fun': Con130}
C_S42 = {'type': 'ineq', 'fun': Con131}
C_S43 = {'type': 'ineq', 'fun': Con132}
C_S44 = {'type': 'ineq', 'fun': Con133}
C_S45 = {'type': 'ineq', 'fun': Con134}
C_S46 = {'type': 'eq', 'fun': Con200}
Cons = [
C_D11, C_D12, C_D13, C_D14, C_D15, C_D16, C_D17, C_D18, C_D19,
C_D21, C_D22,
C_D31, C_D32, C_D33, C_D34, C_D35, C_D36, C_D37, C_D38, C_D39,
C_D41, C_D42, C_D43, C_D44, C_D45,
C_S11, C_S12, C_S13, C_S14, C_S15, C_S16, C_S17, C_S18, C_S19,
C_S21,
C_S31, C_S32, C_S33, C_S34, C_S35, C_S36, C_S37, C_S38, C_S39,
C_S41, C_S42, C_S43, C_S44, C_S45, C_S46,
]
solution = minimize(objective, SV, method='SLSQP',
bounds=bnds, constraints=Cons)
dsT, ssT, Qc0, Qp0 = post_process(
solution,
K=K, P0=P0, dc=dc, sc=sc, ds=ds, ss=ss, Q0=Q0, ED=ED, ES=ES, Qc0=Qc0, Qp0=Qp0,
)
dump(
P0=P0, dc=dc, sc=sc, ds=ds, ss=ss,
Q0=Q0, Qc0=Qc0, Qp0=Qp0, dsT=dsT, ssT=ssT,
)
regression_test(
P0=P0, dc=dc, sc=sc, ds=ds, ss=ss,
Q0=Q0, Qc0=Qc0, Qp0=Qp0, dsT=dsT, ssT=ssT,
)
def post_process(
solution: OptimizeResult, K: int,
P0: np.ndarray, dc: np.ndarray, sc: np.ndarray, ds: np.ndarray, ss: np.ndarray,
Q0: np.ndarray, ED: np.ndarray, ES: np.ndarray, Qc0: np.ndarray, Qp0: np.ndarray,
):
Xx = solution.x
# Return Model Results to Variables
for k in range(K):
if k == 0:
P0[k] = Xx[k]
dc[k] = Xx[1 + 5*k]
sc[k] = Xx[2 + 5*k]
ds[k] = Xx[3 + 5*k]
ss[k] = Xx[4 + 5*k]
Q0[k] = Xx[5 + 5*k]
ED[K] = -1.
dsT = 1/ED * P0/Q0
ds[K] = (P0[K] - dc[K]) / Q0[K]
ES[K] = 1.
ssT = 1/ES * P0/Q0
ss[K] = (P0[K] - sc[K]) / Q0[K]
Qc0 = Qc0[:K]
Qp0 = Qp0[:K]
return dsT, ssT, Qc0, Qp0
def dump(
P0: float, dc: float, sc: float, ds: float, ss: float,
Q0: float, Qc0: float, Qp0: float, dsT: float, ssT: float,
) -> None:
print("Initial Prices")
print(P0)
print()
print("Initial Demand Choke Price")
print(dc)
print("Initial Supply Choke Prices")
print(sc)
print()
print("Initial Slope of the Demand Curve")
print(ds)
print()
print("Initial Slope of the Supply Curve")
print(ss)
print()
print("Initial Quantity")
print(Q0)
print()
print("Initial Quantity Consumed")
print(Qc0)
print()
print("Initial Quantity Produced")
print(Qp0)
print()
print()
print("Initial Slope of the Demand Curve")
print(ds)
print("Test Slope of the Demand Curve")
print(dsT)
print()
print()
print("Initial Slope of the Supply Curve")
print(ss)
print("Test Slope of the Supply Curve")
print(ssT)
print()
def regression_test(
P0: float, dc: float, sc: float, ds: float, ss: float,
Q0: float, Qc0: float, Qp0: float, dsT: float, ssT: float,
) -> None:
assert np.allclose( P0, ( 73.06327500, 51.69400000, 96.94580000, 61.74040000, 73.06327500))
assert np.allclose( dc, ( 131.49283939, 127.71458824, 286.10833659, 116.46024401, 286.10833659))
assert np.allclose( sc, ( 45.48494444, 5.87993628, 39.54745779, 25.18598911, 5.87993628))
assert np.allclose( ds, ( -0.00048987, -0.00577752, -0.01807573, -0.00021679, -0.00056463))
assert np.allclose( ss, ( 0.00035924, 0.00348184, 0.00548479, 0.00014482, 0.00017805))
assert np.allclose(dsT, ( -0.00043716, -0.00577752, -0.01807573, -0.00021679, -0.00019364))
assert np.allclose(ssT, ( 0.00032058, 0.00348184, 0.00548479, 0.00014482, 0.00019364))
assert np.allclose( Q0, (101291.71210000, 13158.00000000, 10465.00000000, 252406.28780000, 377320.9999))
assert np.allclose(Qc0, ( 99861.98900000, 13164.00000000, 10490.00000000, 84311.51100000))
assert np.allclose(Qp0, (102721.43520000, 13152.00000000, 10440.00000000, 420501.06460000))
if __name__ == '__main__':
main()