from pyomo.environ import *

"""
Explanation: Pyomo - Getting started
Pyomo installation: see http://www.pyomo.org/installation
pip install pyomo
End of explanation
"""


model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)

model.obj = Objective(expr = sum(c[n] * model.x[n] for n in N) )

m = 0
model.const_0 = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m] )
m = 1
model.const_1 = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m] )
m = 2
model.const_2 = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m] )

model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail

"""
Explanation: Definition of the problem
$$
\begin{align}
    \min_{x_0,x_1} & \quad -x_0 + 4 x_1 \
    \text{s.t.}    & \quad -3 x_0 + x_1 \leq 6 \
                   & \quad -x_0 - 2 x_1 \geq -4 \
                   & \quad  x_1 \geq -3
\end{align}
$$
Solution
Version 1: naive implementation
End of explanation
"""


model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)

model.obj = Objective(expr = sum(c[n] * model.x[n] for n in N) )

for m in M:
    name = "const_{}".format(m)
    val = Constraint(expr = sum(A[m][n] * model.x[n] for n in N) <= b[m])
    model.add_component(name=name, val=val)

model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail

"""
Explanation: Version 2: without construction rules
End of explanation
"""


model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)

model.obj = Objective(expr = sum(c[n] * model.x[n] for n in N) )

def constraint_fn(model, m):
    return sum(A[m][n] * model.x[n] for n in N) <= b[m]

model.constraint = Constraint(M, rule=constraint_fn)

model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail

"""
Explanation: Version 3: with construction rules
End of explanation
"""


model = ConcreteModel(name="Getting started")

c = [-1, 4]
b = [ 6, 4, 3]
A = [[-3,  1],
     [ 1,  2],
     [ 0, -1]]

N = list(range(len(c)))    # num decision variables
M = list(range(len(A)))    # num constraints

model.x = Var(N)


def objective_fn(model):
    return sum(c[n] * model.x[n] for n in N)

model.obj = Objective(rule=objective_fn)


def constraint_fn(model, m):
    return sum(A[m][n] * model.x[n] for n in N) <= b[m]

model.constraint = Constraint(M, rule=constraint_fn)


model.pprint()

# @tail:
print()
print("-" * 60)
print()

opt = SolverFactory('glpk')

results = opt.solve(model)

model.display()

print()
print("Optimal solution: ", [value(model.x[n]) for n in N])
print("Gain of the optimal solution: ", value(model.obj))
# @:tail

"""
Explanation: Version 4: construction rules on Objective too
End of explanation
"""


