The title may be poorly worded. What I need is coding that expands or contracts according to the number of "countries" included for analysis. The problem requires a minimum of 2 countries but can include N countries depending upon the problem addressed. Using scipy.optimize I can code for 2 countries or for 5 countries, for example, but I cannot figure out how to write code so that the same set of instructions can accommodate 2 countries in the first pass and then 5 countries in the second pass. If all my equations were linear, I can easily add to or subtract from the X matrix and Y vector to solve np.linalg.inv(X).dot(Y) for any number of countries. I want to do the same using scipy.optimize because my equations are non-linear. Note that I am not struck on scipy.optimize. I will use any module as long as it solves the problem and allows for the problem to be expanded or contracted simply by being given data for N countries.
The following coding works just fine. But say I want Trade=K=6 instead of 4. How do I write this program so that K can be solved for K = 4 or 6 without having to rewrite the program?
Also, I am an old fart who learned programming by taking a Fortran class. I think that my python coding looks chunky. I will accept any suggestions as to how to this code more elegant.
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
I'm going to stop short of a full demonstration and just write a list of things that need to change. Do the heavy lifting, and if you have something that doesn't work post a new question (and link it here); if you have something that does work post on codereview.stackexchange.com instead.
float(81.8)is just81.8. Python has a weak type system and adding casts will not help.Why are you adding one element to half of your vectors? If it's truly a "starting value" then it doesn't belong with your other data. It should be used once, for the
x0parameter tominimize, and then thrown away.Get rid of all of your constant-index assignments for input data at the top. Write this CSV, or one like it:
then load it with
pd.read_csv.Your objective and every one of your constraints is wrong. Pay special attention to the signature for
objective(x). That(x)is not just decoration - those are your changing variables, and you're currently ignoring them. You need to unpack your actual decision variables fromx, notXx. To enforce this and discourage accidental closure binding, I strongly encourage you to move all of your code to functions with no global state.Replace
np.delete(P0, [K])with a sliceP0[:-1]. Better yet, those two probably shouldn't be mixed in the first place.For any constraints that support it, you should drop function form and instead use a
LinearConstraint. The only constraints that should be functions are non-linear ones; they should useNonlinearConstraint, and they should be given sibling jacobians rather than letting them be inferred.Don't write line continuation
\on the inside of a parenthesized expression; that isn't necessary.Add regression tests.
An edit doing a small subset of the above can look like