# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.13.6
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# %% [markdown]
# ---
# # <center> <font color=purple>[Pyomo](http://www.pyomo.org/)</font>  </center>
#
# ### <center>A Python-based, open-source optimization modeling language</center>
#
# ---
#
# + [Pyomo documentation](https://pyomo.readthedocs.io/en/stable/)
# + Pyomo is a modeling language, not a solver. It interfaces to solvers such as [glpk](https://www.gnu.org/software/glpk/)
#
# ## Standard import

# %%
import pyomo.environ as pyo
import pandas as pd
import numpy as np

# %% [markdown]
# # [The Diet Problem](https://notebook.community/Pyomo/PyomoGallery/diet/DietProblem)
#
# The goal of the Diet Problem is to select foods that satisfy daily nutritional requirements at minimum cost.
#
#
#
# ### Sets
#  
# | Symbol | Definition |
# | :----- | ---------------: | 
# |  $F$   | set of foods     |
# |  $N$   | set of nutrients |
#
# ### Parameters
#
# | Symbol | Definition                                              | Size               |
# | :----- | :------------------------------------------------------ | :---------------- | 
# | $$c_i$$  | cost per serving of food $i$ | $$\forall i \in F$$ |
# | $$a_{ij}$$   |  amount of nutrient $j$ in food $i$ | $$\forall i \in F, \forall j \in N$$ |
# | $$Nmin_j$$  | minimum level of nutrient $j$  |  $$\forall j \in N$$   |
# | $$Nmax_j$$  | maximum level of nutrient $j$  | $$\forall j \in N$$ |
# | $$V_i$$  | the volume per serving of food $i$  | $$\forall i \in F$$   |
# | $$Vmax$$  |  maximum volume of food consumed |  |
#
#  
# ### Variables
#
# | Symbol | Definition |
# | :----- | ---------------: | 
# |$$x_i$$ | number of servings of food $i$ to consume |
#
# ### Objective
#
# + Minimize the total cost of the food  
#  $$\min \sum_{i \in F} c_i x_i$$
#
# ### Constraints
#
# + Limit nutrient consumption for each nutrient $j \in N$.  
#  $$Nmin_j \leq \sum_{i \in F} a_{ij} x_i \leq Nmax_j\hspace{3em}\forall j \in N$$
#
# + Limit the volume of food consumed  
#  $$\sum_{i \in F} V_i x_i \leq Vmax$$
#  
# + Consumption lower bound  
#  $$x_i \geq 0\hspace{3em}\forall i \in F$$

# %% [markdown]
# ## Pyomo Modeling Components
#
# + Sets
# + Parameters
# + Variables
# + Objectives
# + Constraints
# + Expressions
# + Suffixes
#
#
#

# %% [markdown]
# ### 0. Load the data

# %%
Vmax = 75.0

# %%
Bounds = pd.DataFrame([['Cal',2000, np.inf],
        ['Carbo'   ,  350  ,  375],
        ['Protein' ,   55 , np.inf],
        ['VitA'    ,  100 , np.inf],
        ['VitC'    ,  100 , np.inf],
        ['Calc'    ,  100  , np.inf],
        ['Iron'    ,  100 , np.inf]], 
       columns=['nutrient','min','max'] )
Bounds = Bounds.set_index('nutrient')     
Bounds

# %%
Foods = pd.DataFrame([[  "Cheeseburger",1.84 ,  4.0,510,     34,     28 ,    15 ,     6 ,   30 ,   20  ],
        [  "Ham Sandwich", 2.19 ,  7.5 ,370,     35,     24 ,    15 ,    10 ,   20 ,   20 ],
        [  "Hamburger",1.84 ,  3.5 ,500,     42,     25 ,     6 ,     2 ,   25 ,   20 ],
        [  "Fish Sandwich" , 1.44 ,  5.0 ,370,     38,     14 ,     2 ,     0 ,   15 ,   10 ],
        [  "Chicken Sandwich", 2.29,7.3,400,     42,     31 ,     8 ,    15 ,   15 ,    8],
        [  "Fries"  , .77 ,  2.6 ,220,     26,      3 ,     0 ,    15 ,    0 ,    2 ],
        [  "Sausage Biscuit", 1.29 ,  4.1 ,345,     27,     15 ,     4 ,     0 ,   20 ,   15 ],
        [  "Lowfat Milk" , .60,8.0 ,110,     12,      9 ,    10 ,     4 ,   30 ,    0],
        [  "Orange Juice", .72 , 12.0 ,80 ,    20 ,     1  ,    2  ,  120  ,   2  ,   2 ]],
                 columns = ['food','c','V','Cal' , 'Carbo','Protein' ,  'VitA'  , 'VitC' , 'Calc' , 'Iron'])

Foods.set_index('food' ,inplace=True)
Foods

# %% [markdown]
# ### 1. Create a model object.

# %%
model = pyo.ConcreteModel()

# %% [markdown]
# ### 2. Define sets

# %%
# foods set
model.F = pyo.Set(initialize = Foods.index)

# nutrients set
model.N = pyo.Set(initialize = Bounds.index)    

# %% [markdown]
# ### 3. Define parameters

# %%
model.c       = pyo.Param(model.F,initialize=Foods['c'].to_dict())
model.V       = pyo.Param(model.F,initialize=Foods['V'].to_dict())
model.Cal     = pyo.Param(model.F,initialize=Foods['Cal'].to_dict())
model.Carbo   = pyo.Param(model.F,initialize=Foods['Carbo'].to_dict())
model.Protein = pyo.Param(model.F,initialize=Foods['Protein'].to_dict())
model.VitA    = pyo.Param(model.F,initialize=Foods['VitA'].to_dict())
model.VitC    = pyo.Param(model.F,initialize=Foods['VitC'].to_dict())
model.Calc    = pyo.Param(model.F,initialize=Foods['Calc'].to_dict())
model.Iron    = pyo.Param(model.F,initialize=Foods['Iron'].to_dict())

# %% [markdown]
# ### 4. Define the decision variables

# %%
model.x = pyo.Var(model.F, within=pyo.NonNegativeIntegers)

# %% [markdown]
# ### 5. Define the objective function

# %%
obj = sum(model.c[food]*model.x[food] for food in model.F)
model.cost = pyo.Objective(rule=obj, sense = pyo.minimize)


# %% [markdown]
# ### 6. Define the constraints

# %% [markdown]
# #### Limit nutrient consumption for each nutrient

# %%
def nutrient_rule(model, nutrient):
    value = sum(Foods.loc[food,nutrient]*model.x[food] for food in model.F)
    return pyo.inequality(Bounds.loc[nutrient,'min'], value, Bounds.loc[nutrient,'max'])

model.bounds = pyo.Constraint(model.N, rule=nutrient_rule)

# %% [markdown]
# #### sdasdfasdf

# %%
expr = sum(Foods.loc[food,'V']*model.x[food] for food in model.F) <= Vmax
model.volume = pyo.Constraint(rule=expr)

# %% [markdown]
# ## 7. Request the dual, ie Lagrange multiplliers (optional)

# %%
model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)

# %% [markdown]
# ## Check the model

# %%
model.pprint()

# %% [markdown]
# ## Solve
#
# Cnnect to a solver with `SolverFactory`. Then run the solver with the `solve` method. This returns an object with information about the solution process. The optimal values of the decision variables are written to the `model` object.

# %%
solver = pyo.SolverFactory('ipopt')  # can also use `glpk`

results = solver.solve(model)

print("The solver returned a status of: " + str(results.solver.status))

# %% [markdown]
# Use the `pyo.value` method to extract the results from the model object.

# %% [markdown]
# ### Extract solution values

# %% [markdown]
# #### Optimal cost

# %%
pyo.value(model.cost)

# %% [markdown]
# #### Print optimal decision variables

# %%
{food:pyo.value(model.x[food]) for food in model.F}

# %% [markdown]
# #### Lagrange multipliers

# %%
for c in model.component_objects(pyo.Constraint):
    print ("\n   Constraint",c)
    for index in c:
        value = round(model.dual[c[index]],5)
        print ("\t{:.5s}\t{:.5f}".format(index if index!=None else "None",value))
