Ejercicio 1:¶

In [1]:
!pip install glpk

Collecting glpk
  Downloading glpk-0.4.8.tar.gz (160 kB)
     -------------------------------------- 160.7/160.7 kB 1.9 MB/s eta 0:00:00
  Installing build dependencies: started
  Installing build dependencies: finished with status 'done'
  Getting requirements to build wheel: started
  Getting requirements to build wheel: finished with status 'done'
  Preparing metadata (pyproject.toml): started
  Preparing metadata (pyproject.toml): finished with status 'done'
Building wheels for collected packages: glpk
  Building wheel for glpk (pyproject.toml): started
  Building wheel for glpk (pyproject.toml): finished with status 'error'
Failed to build glpk

  error: subprocess-exited-with-error
  
  Building wheel for glpk (pyproject.toml) did not run successfully.
  exit code: 1
  
  [5 lines of output]
  running bdist_wheel
  running build
  running build_ext
  building 'glpk' extension
  error: Microsoft Visual C++ 14.0 or greater is required. Get it with "Microsoft C++ Build Tools": https://visualstudio.microsoft.com/visual-cpp-build-tools/
  [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
  ERROR: Failed building wheel for glpk
ERROR: Could not build wheels for glpk, which is required to install pyproject.toml-based projects
In [2]:
!pip install pulp

Collecting pulp
  Using cached PuLP-2.9.0-py3-none-any.whl (17.7 MB)
Installing collected packages: pulp
Successfully installed pulp-2.9.0
In [2]:
import pulp as pp
In [3]:
model = pp.LpProblem(name='medicine-problem', # just the name
                     sense=pp.LpMinimize) # type of problem
In [4]:
# how much gas?
VV = pp.LpVariable(name="VegaVita",  # just the name
                    lowBound=0,  # ensure non-negativity
                    cat='Continuous') # here: you accept decimal values

# how much oil?
HH = pp.LpVariable(name="HappyHealth",
                 lowBound=0,
                 cat='Continuous')
In [5]:
VVCoeff=0.2
HHCoeff=1.5
obj_func = VVCoeff*VV + HHCoeff*HH
In [6]:
C1= pp.LpConstraint(name='Vitamin C Constraint',   # just the name
                    e= 20*VV + 30*HH, rhs=60, # linear combination of constraint and rhs
                    sense=pp.LpConstraintGE) # 'rule' >= 0 (LpConstraintGE)
C2= pp.LpConstraint(name='Calcium Constraint',
                    e= 500*VV + 250*HH, rhs=1000,
                    sense=pp.LpConstraintGE) # 'rule' >= 3000000 (LpConstraintGE)
C3= pp.LpConstraint(name='Iron Constraint',
                    e= 9*VV + 2*HH, rhs=18,
                    sense=pp.LpConstraintGE)
C4= pp.LpConstraint(name='Niacin Constraint',
                    e= 2*VV + 10*HH, rhs=20,
                    sense=pp.LpConstraintGE)
C5= pp.LpConstraint(name='Mg Constraint',
                    e= 60*VV + 90*HH, rhs=360,
                    sense=pp.LpConstraintGE, )# 'rule' <= 6400000 (LpConstraintLE)
In [7]:
model += obj_func
model += C1
model += C2
model += C3
In [8]:
solver_list = pp.listSolvers()
print(solver_list)

['GLPK_CMD', 'PYGLPK', 'CPLEX_CMD', 'CPLEX_PY', 'GUROBI', 'GUROBI_CMD', 'MOSEK', 'XPRESS', 'XPRESS', 'XPRESS_PY', 'PULP_CBC_CMD', 'COIN_CMD', 'COINMP_DLL', 'CHOCO_CMD', 'MIPCL_CMD', 'SCIP_CMD', 'FSCIP_CMD', 'SCIP_PY', 'HiGHS', 'HiGHS_CMD', 'COPT', 'COPT_DLL', 'COPT_CMD']
In [9]:
solverToUse= pp.COIN_CMD(msg=False)
model.solve();
In [10]:
import pandas as pd

Results={"Model Status":pp.LpStatus[model.status]}
Results.update({"Optimal Solution":pp.value(model.objective)})
Results.update({v.name: v.varValue for v in model.variables()})
Results
Out[10]:
{'Model Status': 'Optimal',
 'Optimal Solution': 0.6000000000000001,
 'HappyHealth': 0.0,
 'VegaVita': 3.0}

Ejercicio 2¶

Exercise: Choosing a country for a Master Program¶

  • Make a group of 4 people from this course.
  • If you have the criteria: cost of living, language difficulty, possibilities to get a job in that country after studies are finished.
  • If you have the alternatives: Brazil, Spain, USA, Germany.
  • Create an AHP model and get the ranking.

We choose¶

  • Master Program: A Master of Public Policy (MPP)
  • Criteria: (1) cost of living, (2) language difficulty, (3) possibilities to get a job in that country after studies are finished.
  • Alternatives: UK, Germany, Spain and Chile.
  • AHP model: https://docs.google.com/spreadsheets/d/1fcnfcmDSSisiYiu-FpnYT5cKBQYqrjnE1DkigsHcM8E/edit?gid=659603613#gid=659603613

AHP Model¶

  1. Prepare data for comparision
In [58]:
#%%html

#<iframe src="https://docs.google.com/spreadsheets/d/1fcnfcmDSSisiYiu-FpnYT5cKBQYqrjnE1DkigsHcM8E/edit?gid=659603613#gid=659603613" width="600" height="300" ></iframe>
In [57]:
%%html

<iframe src="https://docs.google.com/spreadsheets/d/1hwWuzUFCzs_YFoV27HWYzSMLgkbXCpIc/edit?gid=113328570#gid=113328570" width="600" height="300" ></iframe>
  1. Get the data (Excel) and Open each sheet:
In [59]:
#linkGoogle='https://docs.google.com/spreadsheets/d/e/2PACX-1vRGzlHwDqOoeb337eloDICBNAUUVdKz47f2drs80oAB8Xhba900zF4dNFBhIOjOyg/pub?output=xlsx'# the link to the data
In [67]:
import pandas as pd

# Enlace directo de Google Sheets en formato .xlsx
linkGoogle = 'https://docs.google.com/spreadsheets/d/1hwWuzUFCzs_YFoV27HWYzSMLgkbXCpIc/export?format=xlsx&id=1hwWuzUFCzs_YFoV27HWYzSMLgkbXCpIc'

# Leer las diferentes hojas
pairwise_cost = pd.read_excel(linkGoogle, sheet_name='cost of living', index_col=0, engine='openpyxl')
pairwise_language = pd.read_excel(linkGoogle, sheet_name='languaje difficulty', index_col=0, engine='openpyxl')
pairwise_job = pd.read_excel(linkGoogle, sheet_name='possibilities to get a job', index_col=0, engine='openpyxl')
pairwise_criterio = pd.read_excel(linkGoogle,sheet_name='criterios', index_col=0, engine='openpyxl')
In [68]:
#check
pairwise_criterio
Out[68]:
cost language job
cost 1.0 5 0.111111
language 0.2 1 0.142857
job 9.0 7 1.000000
In [69]:
import networkx as nx

cost = nx.from_pandas_adjacency(pairwise_cost,create_using=nx.MultiDiGraph())

# pairwise
cost.edges(data=True)
Out[69]:
OutMultiEdgeDataView([('uk', 'uk', {'weight': 1.0}), ('uk', 'germany', {'weight': 0.3333}), ('uk', 'spain', {'weight': 0.5}), ('uk', 'chile', {'weight': 0.11111}), ('germany', 'germany', {'weight': 1.0}), ('germany', 'uk', {'weight': 3.0}), ('germany', 'spain', {'weight': 0.5}), ('germany', 'chile', {'weight': 0.2}), ('spain', 'spain', {'weight': 1.0}), ('spain', 'uk', {'weight': 2.0}), ('spain', 'germany', {'weight': 2.0}), ('spain', 'chile', {'weight': 0.5}), ('chile', 'chile', {'weight': 1.0}), ('chile', 'uk', {'weight': 9.0}), ('chile', 'germany', {'weight': 5.0}), ('chile', 'spain', {'weight': 2.0})])
In [70]:
cost_comparisons ={(e[0],e[1]):e[2]['weight'] for e in cost.edges(data=True) if e[0]!= e[1]}
cost_comparisons
Out[70]:
{('uk', 'germany'): 0.3333,
 ('uk', 'spain'): 0.5,
 ('uk', 'chile'): 0.11111,
 ('germany', 'uk'): 3.0,
 ('germany', 'spain'): 0.5,
 ('germany', 'chile'): 0.2,
 ('spain', 'uk'): 2.0,
 ('spain', 'germany'): 2.0,
 ('spain', 'chile'): 0.5,
 ('chile', 'uk'): 9.0,
 ('chile', 'germany'): 5.0,
 ('chile', 'spain'): 2.0}
In [71]:
language = nx.from_pandas_adjacency(pairwise_languaje,create_using=nx.MultiDiGraph())
language_comparisons={(e[0],e[1]):e[2]['weight'] for e in language.edges(data=True) if e[0]!= e[1]}

job = nx.from_pandas_adjacency(pairwise_job,create_using=nx.MultiDiGraph())
job_comparisons={(e[0],e[1]):e[2]['weight'] for e in job.edges(data=True) if e[0]!= e[1]}
In [72]:
[cost_comparisons, language_comparisons,job_comparisons]
Out[72]:
[{('uk', 'germany'): 0.3333,
  ('uk', 'spain'): 0.5,
  ('uk', 'chile'): 0.11111,
  ('germany', 'uk'): 3.0,
  ('germany', 'spain'): 0.5,
  ('germany', 'chile'): 0.2,
  ('spain', 'uk'): 2.0,
  ('spain', 'germany'): 2.0,
  ('spain', 'chile'): 0.5,
  ('chile', 'uk'): 9.0,
  ('chile', 'germany'): 5.0,
  ('chile', 'spain'): 2.0},
 {('uk', 'germany'): 5.0,
  ('uk', 'spain'): 0.1111111111,
  ('uk', 'chile'): 0.1111111111,
  ('germany', 'uk'): 0.2,
  ('germany', 'spain'): 0.1111111111,
  ('germany', 'chile'): 0.1111111111,
  ('spain', 'uk'): 9.0,
  ('spain', 'germany'): 9.0,
  ('spain', 'chile'): 1.0,
  ('chile', 'uk'): 9.0,
  ('chile', 'germany'): 9.0,
  ('chile', 'spain'): 1.0},
 {('uk', 'germany'): 0.33,
  ('uk', 'spain'): 3.0,
  ('uk', 'chile'): 7.0,
  ('germany', 'uk'): 3.0,
  ('germany', 'spain'): 5.0,
  ('germany', 'chile'): 9.0,
  ('spain', 'uk'): 3.0,
  ('spain', 'germany'): 0.2,
  ('spain', 'chile'): 5.0,
  ('chile', 'uk'): 0.14285714285714285,
  ('chile', 'germany'): 0.11,
  ('chile', 'spain'): 0.2}]
In [73]:
G_CRIT = nx.from_pandas_adjacency(pairwise_criterio,create_using=nx.MultiDiGraph())
criteria_comparisons ={(e[0],e[1]):e[2]['weight'] for e in G_CRIT.edges(data=True) if e[0]!= e[1]}
criteria_comparisons
Out[73]:
{('cost', 'language'): 5.0,
 ('cost', 'job'): 0.1111111111111111,
 ('language', 'cost'): 0.2,
 ('language', 'job'): 0.14285714285714285,
 ('job', 'cost'): 9.0,
 ('job', 'language'): 7.0}
In [74]:
!pip install ahpy

Requirement already satisfied: ahpy in c:\users\michael encalada\anaconda3\lib\site-packages (2.0)
Requirement already satisfied: numpy in c:\users\michael encalada\anaconda3\lib\site-packages (from ahpy) (1.24.4)
Requirement already satisfied: scipy in c:\users\michael encalada\anaconda3\lib\site-packages (from ahpy) (1.9.1)
In [75]:
import ahpy
In [76]:
cost_living = ahpy.Compare('cost', cost_comparisons, precision=3, random_index='saaty')
language_difficulty = ahpy.Compare('language', language_comparisons, precision=3, random_index='saaty')
job_opportunities = ahpy.Compare('job', job_comparisons, precision=3, random_index='saaty')
criteria = ahpy.Compare('criteria', criteria_comparisons, precision=3, random_index='saaty')
In [77]:
criteria.add_children([cost_living, language_difficulty, job_opportunities])
In [78]:
print(criteria.target_weights)

{'germany': 0.473, 'spain': 0.24, 'chile': 0.148, 'uk': 0.138}
In [79]:
[(val.name,val.consistency_ratio) for val in [cost_living, language_difficulty, job_opportunities, criteria]]
Out[79]:
[('cost', 0.052), ('language', 0.127), ('job', 0.155), ('criteria', 0.382)]

Ejercicio 3¶

In [80]:
import os

print(os.getcwd())  # Esto muestra el directorio actual
print(os.listdir()) # Esto lista todos los archivos en el directorio actual

C:\Users\Michael Encalada\Documents\GitHub\Magallanes\week12_optimization
['.git', '.gitignore', '.ipynb_checkpoints', 'Ejercicios.html', 'Ejercicios.ipynb', 'LICENSE', 'municipalidades.xlsx', 'Optimization tarea.xlsx', 'README.md']
In [81]:
import pandas as pd

df = pd.read_excel(r"C:\Users\Michael Encalada\Documents\GitHub\Magallanes\week12_optimization\municipalidades.xlsx")  # Usa r"" para indicar una cadena de texto "raw" y evitar problemas de escape en Windows
In [82]:
df.head()  # Muestra las primeras filas del DataFrame
Out[82]:
muni Ingreso per capita % Presupuesto ejecutado trabajadores n° beneficiarios beneficiarios% centros años implementación
0 ANCON 207.857758 68.46 215 50 0.030395 1 10
1 ATE 289.322960 77.10 1627 420 0.021511 12 5
2 BARRANCO 292.157317 75.90 521 45 0.007888 1 0
3 CHACLACAYO 213.441099 80.92 189 180 0.052988 4 11
4 CHORRILLOS 264.977511 83.82 345 361 0.021194 3 13
In [83]:
# Calcular los ratios con los nombres de columna adecuados
df['rate1'] = df['n° beneficiarios'] / df['trabajadores']  # Ratio de beneficiarios por trabajador
df['rate2'] = df['% Presupuesto ejecutado'] / df['Ingreso per capita']  # Ratio de presupuesto ejecutado relativo al ingreso per cápita

# Visualizar los resultados
print(df[['muni', 'rate1', 'rate2']])

                         muni     rate1     rate2
0                       ANCON  0.232558  0.329360
1                         ATE  0.258144  0.266484
2                    BARRANCO  0.086372  0.259792
3                  CHACLACAYO  0.952381  0.379121
4                  CHORRILLOS  1.046377  0.316329
5   LOS OLIVOS (LAS PALMERAS)  0.392377  0.308109
6        LURIGANCHO (CHOSICA)  0.068729  0.292089
7               PUENTE PIEDRA  0.856531  0.349001
8                    SAN LUIS  1.157407  0.334111
9                  SANTA ROSA  0.808081  0.275528
10                  SURQUILLO  1.287554  0.356133
In [84]:
!pip install altair

Requirement already satisfied: altair in c:\users\michael encalada\anaconda3\lib\site-packages (5.4.1)
Requirement already satisfied: jsonschema>=3.0 in c:\users\michael encalada\anaconda3\lib\site-packages (from altair) (4.16.0)
Requirement already satisfied: narwhals>=1.5.2 in c:\users\michael encalada\anaconda3\lib\site-packages (from altair) (1.13.3)
Requirement already satisfied: packaging in c:\users\michael encalada\anaconda3\lib\site-packages (from altair) (21.3)
Requirement already satisfied: typing-extensions>=4.10.0 in c:\users\michael encalada\anaconda3\lib\site-packages (from altair) (4.12.2)
Requirement already satisfied: jinja2 in c:\users\michael encalada\anaconda3\lib\site-packages (from altair) (3.1.2)
Requirement already satisfied: attrs>=17.4.0 in c:\users\michael encalada\anaconda3\lib\site-packages (from jsonschema>=3.0->altair) (21.4.0)
Requirement already satisfied: pyrsistent!=0.17.0,!=0.17.1,!=0.17.2,>=0.14.0 in c:\users\michael encalada\anaconda3\lib\site-packages (from jsonschema>=3.0->altair) (0.18.0)
Requirement already satisfied: MarkupSafe>=2.0 in c:\users\michael encalada\anaconda3\lib\site-packages (from jinja2->altair) (3.0.2)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in c:\users\michael encalada\anaconda3\lib\site-packages (from packaging->altair) (3.0.9)
In [85]:
import altair as alt

points = alt.Chart(df).mark_point().encode(
    x=alt.X('rate1:Q', title='Beneficiarios vs Trabajadores'),  # Cambia el título de este eje
    y=alt.Y('rate2:Q', title='Presupuesto ejecutado vs Ingresos, per cápita')   # Cambia el título de este eje
)

text = points.mark_text(
    align='right',
    baseline='middle',
    dx=-7
).encode(
    text='muni'
).interactive()

points + text
Out[85]:
In [86]:
df[['muni','rate1','rate2']].sort_values(by='rate1',ascending=False).head()
Out[86]:
muni rate1 rate2
10 SURQUILLO 1.287554 0.356133
8 SAN LUIS 1.157407 0.334111
4 CHORRILLOS 1.046377 0.316329
3 CHACLACAYO 0.952381 0.379121
7 PUENTE PIEDRA 0.856531 0.349001
In [87]:
df[['muni','rate1','rate2']].sort_values(by='rate2',ascending=False).head()
Out[87]:
muni rate1 rate2
3 CHACLACAYO 0.952381 0.379121
10 SURQUILLO 1.287554 0.356133
7 PUENTE PIEDRA 0.856531 0.349001
8 SAN LUIS 1.157407 0.334111
0 ANCON 0.232558 0.329360
In [88]:
Best_beneficiarios=df.rate1.idxmax()
Best_dinero=df.rate2.idxmax()

frontier1=df.loc[Best_beneficiarios,['rate1','rate2']].to_list()
frontier2=df.loc[Best_dinero,['rate1','rate2']].to_list()

#parallels
frontier1v=[frontier1[0],0]
frontier2h=[0,frontier2[1]]

#then
envelope=pd.DataFrame([frontier2h,frontier2,frontier1,frontier1v],columns=['x','y'])
envelope
Out[88]:
x y
0 0.000000 0.379121
1 0.952381 0.379121
2 1.287554 0.356133
3 1.287554 0.000000
In [89]:
points + text + alt.Chart(envelope).mark_line(color='red').encode(
    x='x',
    y='y',
)
Out[89]:

Interpretación La línea roja representa un límite de eficiencia o frontera de producción, que suele utilizarse en análisis de eficiencia como el Análisis Envolvente de Datos (DEA). Beneficiarios trabajadores: Los valores en el eje X muestran cómo los municipios varían en términos de la relación entre beneficiarios y trabajadores. Los municipios más a la derecha tienen una mayor relación "beneficiarios/trabajadores", lo que sugiere que están atendiendo a más beneficiarios por trabajador. Presupuesto vs ingresos per cápita Los valores en el eje Y representan la relación entre ingresos y presupuesto ejecutado. Los municipios que están más arriba en el eje Y pueden estar mostrando una mayor relación de presupuesto ejecutado per capita respecto a sus ingresos, lo que podría indicar una mejor gestión de sus recursos financieros en comparación con los que están más abajo.

In [90]:
!pip install Pyfrontier

Requirement already satisfied: Pyfrontier in c:\users\michael encalada\anaconda3\lib\site-packages (1.0.2)
Requirement already satisfied: pulp>=2.6.0 in c:\users\michael encalada\anaconda3\lib\site-packages (from Pyfrontier) (2.9.0)
Requirement already satisfied: numpy>=1.21 in c:\users\michael encalada\anaconda3\lib\site-packages (from Pyfrontier) (1.24.4)
In [91]:
Input=df.columns[1:2].to_list()
Output=df.columns[3:7].to_list()
In [92]:
from Pyfrontier.frontier_model import EnvelopDEA

dea_air_vrs_in = EnvelopDEA("VRS", "in")
dea_air_vrs_in.fit(
    inputs=df[Input].to_numpy(),
    outputs=df[Output].to_numpy()
)
In [93]:
df['vrs_in']=[r.score for r in dea_air_vrs_in.result]
df.set_index(df.muni,inplace=True)
df['vrs_in']
Out[93]:
muni
ANCON                        1.000000
ATE                          1.000000
BARRANCO                     0.753817
CHACLACAYO                   1.000000
CHORRILLOS                   0.817271
LOS OLIVOS (LAS PALMERAS)    1.000000
LURIGANCHO (CHOSICA)         0.936022
PUENTE PIEDRA                1.000000
SAN LUIS                     0.881591
SANTA ROSA                   1.000000
SURQUILLO                    1.000000
Name: vrs_in, dtype: float64
In [94]:
!pip install py4etrics

Requirement already satisfied: py4etrics in c:\users\michael encalada\anaconda3\lib\site-packages (0.1.9)
In [95]:
import numpy as np # linear algebra
from py4etrics import tobit
# import statsmodels.api as sm

df['censored'] =np.where(df['vrs_in']==1, 1, 0)
cens = df['censored']
endog = df.loc[:,'vrs_in']
exog = df.loc[:,'trabajadores':'n° beneficiarios']

tobit_res = tobit.Tobit(endog, exog, cens,right=1).fit()
tobit_res.summary()

Optimization terminated successfully.
         Current function value: 0.733342
         Iterations: 83
         Function evaluations: 151
Out[95]:
Tobit Regression Results
Dep. Variable: vrs_in Pseudo R-squ: -1.602
Method: Maximum Likelihood Log-Likelihood: -8.1
No. Observations: 11 LL-Null: -3.1
No. Uncensored Obs: 4 LL-Ratio: -9.9
No. Left-censored Obs: 0 LLR p-value: 1.000
No. Right-censored Obs: 7 AIC: 18.1
Df Residuals: 9 BIC: 18.5
Df Model: 1 Covariance Type: nonrobust
coef std err z P>|z| [0.025 0.975]
trabajadores 0.0013 0.001 1.769 0.077 -0.000 0.003
n° beneficiarios 0.0019 0.001 1.447 0.148 -0.001 0.005
Log(Sigma) -0.2939 0.394 -0.747 0.455 -1.066 0.478

Resultados: En atención al adulto mayor, y los beneficiarios de servicios ofrecidos por los municipios se observa que el municipio "SURQUILLO" podrían estar considerados como más eficientes, ya que optimizan su relación entre beneficiarios y trabajadores mientras mantienen una buena proporción de ingresos sobre el presupuesto.

Municipios alejados de la línea roja, como "PUENTE PIEDRA" y "CHACLACAYO", podrían tener margen para mejorar su eficiencia en alguna de las dimensiones que estás evaluando.

FIN¶

In [ ]: