連立方程式体系を最適化問題に変換する(1)

前回の記事の続き。


import pyomo.environ as pyo


model = pyo.ConcreteModel(name="utility maximization",
                          doc="8 variables, 8 constraints")

model.x1a = pyo.Var(domain=pyo.NonNegativeReals)
model.x2a = pyo.Var(domain=pyo.NonNegativeReals)
model.x1b = pyo.Var(domain=pyo.NonNegativeReals)
model.x2b = pyo.Var(domain=pyo.NonNegativeReals)
model.x1 = pyo.Var(domain=pyo.NonNegativeReals)
model.x2 = pyo.Var(domain=pyo.NonNegativeReals)
model.p1 = pyo.Var(domain=pyo.NonNegativeReals)
model.p2 = pyo.Var(domain=pyo.NonNegativeReals)


model.OBJ = pyo.Objective(expr = model.x1**2 * model.x2,
                          sense = pyo.maximize)

model.Constraint1 = pyo.Constraint(expr = model.x1a == 140+25*model.p1**(-2)*model.p2**2)
model.Constraint2 = pyo.Constraint(expr = model.x2a == 70*model.p1*model.p2**(-1)+12.5*model.p1**(-1)*model.p2)
model.Constraint3 = pyo.Constraint(expr = model.x1b == 190+12.5*model.p1**(-2)*model.p2**2)
model.Constraint4 = pyo.Constraint(expr = model.x2b == 380*model.p1*model.p2**(-1)+25*model.p1**(-1)*model.p2)
model.Constraint5 = pyo.Constraint(expr = model.x2 == 17.32*model.x1**0.5)
model.Constraint6 = pyo.Constraint(expr = model.x1 == 75*model.p1**(-2)*model.p2**2)
model.Constraint7 = pyo.Constraint(expr = model.x1a +model.x1b +model.x1==780)
model.Constraint8 = pyo.Constraint(expr = model.x2a +model.x2b== model.x2)



opt = pyo.SolverFactory('ipopt')
res = opt.solve(model)

print(f"評価関数:{model.OBJ()}")
print(f"x1: {model.x1()}")
print(f"x2: {model.x2()}")