Economic equilibria problem with PATH solver#
Description: showcase of the PATH solver via Kestrel interface and Neos server to solve economic model using complementarity conditions from the AMPL book
Tags: ampl-book, PATH, Kestrel, Neos, finance, complementarity-problem
Notebook author: Marcos Dominguez Velad <marcos@ampl.com>
Model author: N/A
# Install dependencies
%pip install -q amplpy pandas
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["gokestrel"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
Context#
The main purpose of this notebook is to show how to use the PATH solver with AMPL.
Currently, one of the ways to access it seamlessly is via the Neos server using the Kestrel interface, which is open-source and available through a Community Edition.
Kestrel setup#
As showed at the beginning of the notebook, these are the steps to get Kestrel (gokestrel module) working:
# Install Python API for AMPL:
$ python -m pip install amplpy --upgrade
# Install AMPL & solver modules:
$ python -m amplpy.modules install gokestrel # install Kestrel
# Activate your license (e.g., free ampl.com/ce or ampl.com/courses licenses):
$ python -m amplpy.modules activate <your-license-uuid>
Or simply, in a Jupyter notebook.
# Install dependencies
%pip install -q amplpy pandas
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["gokestrel"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
The Kestrel interface will use AMPL locally with remote solvers from Neos:
from amplpy import AMPL
ampl = AMPL()
...
ampl.option = {
'solver' : 'kestrel',
'email' : 'your-email-address@xxx.com',
'kestrel_options' : 'solver=PATH, ...',
}
ampl.solve()
Many other solvers are available with Kestrel.
Economic equilibria#
This model is based on the models showed in the Chapter 19 of the AMPL book, related to complementarity problems.
Classic model#
The following model does not use complementarity conditions.
Sets and parameters:
Products
PRODanddemand[i]for each producti.Activities
ACTandcost[j]for each activityj.Amount of product
iproduced by activityj,io[i,j].
Decision variables
Level[j]: levels of production activities (for each activityj).Objective: minimize total production cost,
\[
\sum \limits_{j \in ACT} cost_j \cdot Level_j;
\]
So the AMPL formulation is implemented as:
%%writefile econmin.mod
set PROD; # products
set ACT; # activities
param cost {ACT} > 0; # cost per unit of each activity
param demand {PROD} >= 0; # units of demand for each product
param io {PROD,ACT} >= 0; # units of each product from
# 1 unit of each activity
var Level {j in ACT} >= 0;
minimize Total_Cost: sum {j in ACT} cost[j] * Level[j];
subject to Demand {i in PROD}:
sum {j in ACT} io[i,j] * Level[j] >= demand[i];
Overwriting econmin.mod
import pandas as pd
# ACT cost data
act_cost_data = {
"ACT": ["P1", "P1a", "P2", "P2a", "P2b", "P3", "P3c", "P4"],
"cost": [2450, 1290, 1850, 3700, 2150, 2200, 2370, 2170],
}
act_cost_df = pd.DataFrame(act_cost_data).set_index("ACT")
# PROD demand data
prod_demand_data = {
"PROD": ["AA1", "AC1", "BC1", "BC2", "NA2", "NA3"],
"demand": [70000, 80000, 90000, 70000, 400000, 800000],
}
prod_demand_df = pd.DataFrame(prod_demand_data).set_index("PROD")
# IO matrix data
io_data = {
"AA1": [60, 8, 8, 40, 15, 70, 25, 60],
"AC1": [20, 0, 10, 40, 35, 30, 40, 20],
"BC1": [10, 20, 15, 35, 15, 15, 30, 15],
"BC2": [15, 20, 10, 10, 15, 15, 30, 10],
"NA2": [938, 1180, 945, 278, 1182, 896, 1029, 1397],
"NA3": [295, 770, 440, 430, 315, 400, 370, 450],
}
io_index = ["P1", "P1a", "P2", "P2a", "P2b", "P3", "P3c", "P4"]
io_df = pd.DataFrame(io_data, index=io_index)
ampl.reset()
ampl.read("econmin.mod")
ampl.set_data(prod_demand_df, "PROD")
ampl.set_data(act_cost_df, "ACT")
ampl.param["io"] = io_df.T
ampl.option = {
"solver": "kestrel",
"email": "your-email-address@xxx.com",
"kestrel_options": "solver=PATH",
}
ampl.solve()
ampl.display("Level, Demand.dual")
Connecting to: neos-server.org:3333
Job 17536349 submitted to NEOS, password='oKjEqriZ'
Check the following URL for progress report:
https://neos-server.org/neos/cgi-bin/nph-neos-solver.cgi?admin=results&jobnumber=17536349&pass=oKjEqriZ
Job 17536349 dispatched
password: oKjEqriZ
---------- Begin Solver Output -----------
Condor submit: 'neos.submit'
Condor submit: 'watchdog.submit'
Job submitted to NEOS HTCondor pool.
kestrel_options:solver=path
Supplied solver options are: threads=4
You are using the solver path.
Executing on prod-exec-6.neos-server.org
Path 4.7.03: threads=4
14 side inequalities:
8 simple bounds,
0 range bounds,
6 single inequality constraints, and
0 range constraints.
Path 4.7.03 (Sat May 31 16:42:27 2014)
Written by Todd Munson, Steven Dirkse, and Michael Ferris
INITIAL POINT STATISTICS
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[1])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[2])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[3])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[4])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[5])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[6])
Zero row of order . . . . . . . 0.0000e+00 eqn: (**con_name(bad n)**)
Zero row of order . . . . . . . 0.0000e+00 eqn: (**con_name(bad n)**)
Total zero rows . . . . . . . . 8
Maximum of X. . . . . . . . . . 0.0000e+00 var: (_svar[1])
Maximum of F. . . . . . . . . . 8.0000e+05 eqn: (**con_name(bad n)**)
Maximum of Grad F . . . . . . . 1.3970e+03 eqn: (**con_name(bad n)**)
var: (_svar[8])
INITIAL JACOBIAN NORM STATISTICS
Maximum Row Norm. . . . . . . . 7.8460e+03 eqn: (**con_name(bad n)**)
Minimum Row Norm. . . . . . . . 2.0000e+00 eqn: (**con_name(bad n)**)
Maximum Column Norm . . . . . . 1.9990e+03 var: (_svar[2])
Minimum Column Norm . . . . . . 1.0000e+00 var: (**var_name(bad n)**)
Crash Log
major func diff size residual step prox (label)
0 0 1.8158e+06 0.0e+00 (**con_name(bad n)**)
pn_search terminated: no progress.
Major Iteration Log
major minor func grad residual step type prox inorm (label)
0 0 13 1 1.8158e+06 I 1.0e+01 1.6e+06 (**con_name(bad )
1 1 14 2 2.1557e+06 1.0e+00 SO 4.0e+00 1.9e+06 (**con_name(bad )
2 1 15 3 4.8834e+06 1.0e+00 SO 1.6e+00 4.3e+06 (**con_name(bad )
3 2 16 4 1.8158e+06 1.0e+00 RO 6.4e-01 1.6e+06 (**con_name(bad )
4 2 17 5 1.8158e+06 1.0e+00 RO 2.6e-01 1.6e+06 (**con_name(bad )
5 2 18 6 1.8158e+06 1.0e+00 RO 1.0e-01 1.6e+06 (**con_name(bad )
6 2 19 7 1.8158e+06 1.0e+00 RO 4.1e-02 1.6e+06 (**con_name(bad )
7 2 20 8 1.8158e+06 1.0e+00 RO 1.6e-02 1.6e+06 (**con_name(bad )
8 2 33 9 0.0000e+00 5.4e-05 RG 4.0e+00 0.0e+00 (_scon[1])
FINAL STATISTICS
Inf-Norm of Complementarity . . 0.0000e+00 eqn: (_scon[1])
Inf-Norm of Normal Map. . . . . 0.0000e+00 eqn: (_scon[1])
Inf-Norm of Minimum Map . . . . 0.0000e+00 eqn: (_scon[1])
Inf-Norm of Fischer Function. . 0.0000e+00 eqn: (_scon[1])
Inf-Norm of Grad Fischer Fcn. . 0.0000e+00 eqn: (_scon[1])
Two-Norm of Grad Fischer Fcn. . 0.0000e+00
FINAL POINT STATISTICS
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[1])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[2])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[3])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[4])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[5])
Zero row of order . . . . . . . 0.0000e+00 eqn: (_scon[6])
Zero row of order . . . . . . . 0.0000e+00 eqn: (**con_name(bad n)**)
Zero row of order . . . . . . . 0.0000e+00 eqn: (**con_name(bad n)**)
Total zero rows . . . . . . . . 8
Maximum of X. . . . . . . . . . 2.3776e+05 var: (_svar[2])
Maximum of F. . . . . . . . . . 1.3508e+09 eqn: (**con_name(bad n)**)
Maximum of Grad F . . . . . . . 1.3970e+03 eqn: (**con_name(bad n)**)
var: (_svar[8])
** EXIT - solution found.
Major Iterations. . . . 8
Minor Iterations. . . . 14
Restarts. . . . . . . . 0
Crash Iterations. . . . 0
Gradient Steps. . . . . 1
Function Evaluations. . 33
Gradient Evaluations. . 9
Basis Time. . . . . . . 0.000121
Total Time. . . . . . . 0.000432
Residual. . . . . . . . 0.000000e+00
Path 4.7.03: Solution found.
8 iterations (0 for crash); 14 pivots.
33 function, 9 gradient evaluations.
: Level Demand.dual :=
AA1 . 0
AC1 . 0
BC1 . 0
BC2 . 0
NA2 . 0
NA3 . 0
P1 134795 .
P1a 237764 .
P2 159723 .
P2a 101279 .
P2b 159209 .
P3 149855 .
P3c 156225 .
P4 201806 .
;
assert ampl.solve_result == "solved", ampl.solve_result
Complementarity model#
Consider the new variables Price[i] for each product, and solve the problem to find an equilibrium instead of an optimum solution. The equilibrium is subject to two conditions:
First, for each product, total output must meet demand and price must be nonnegative, and in addition there must be a complementarity between these relationships, where production exceeds demand the price must be zero, or equivalently, if the price is positive, the production must equal the demand.
subject to Pri_Compl {i in PROD}:
Price[i] >= 0 complements
sum {j in ACT} io[i,j] * Level[j] >= demand[i];
Second, for each activity
j, the value of resulting productiisPrice[i]*io[i,j], so activityjwould produce a total value of
\[\sum \limits_{i \in PROD} Price[i] \cdot io[i,j]\]
When equilibrium happens, the previous value must not exceed the activity’s cost per unit cost[j]. Moreover, there is a complementarity between this relationship and the level of activity j, where cost exceeds total value the activity must be zero, or, where the activity cost is positive then the total value must be equal to the cost. This can be expressed as:
subject to Lev_Compl {j in ACT}:
Level[j] >= 0 complements
sum {i in PROD} Price[i] * io[i,j] <= cost[j];
%%writefile econ.mod
set PROD; # products
set ACT; # activities
param cost {ACT} > 0; # cost per unit of each activity
param demand {PROD} >= 0; # units of demand for each product
param io {PROD,ACT} >= 0; # units of each product from
# 1 unit of each activity
var Price {i in PROD};
var Level {j in ACT};
subject to Pri_Compl {i in PROD}:
Price[i] >= 0 complements
sum {j in ACT} io[i,j] * Level[j] >= demand[i];
subject to Lev_Compl {j in ACT}:
Level[j] >= 0 complements
sum {i in PROD} Price[i] * io[i,j] <= cost[j];
Overwriting econ.mod
Remark: the model is square, as there are \(n+m\) variables and \(n+m\) constraints (\(n\) number of products and \(m\) number of activities). Some solvers can take advantage of this to use different solving techniques.
Finally, let’s see that the equilibrium model gives also the optimal solution for the first classic formulation.
ampl.reset()
ampl.read("econ.mod")
ampl.set_data(prod_demand_df, "PROD")
ampl.set_data(act_cost_df, "ACT")
ampl.param["io"] = io_df.T
ampl.option = {
"solver": "kestrel",
"email": "your-email-address@xxx.com",
"kestrel_options": "solver=PATH",
}
ampl.solve()
ampl.display("Level, sum {j in ACT} cost[j] * Level[j]")
Connecting to: neos-server.org:3333
Job 17536350 submitted to NEOS, password='FaxWLjyN'
Check the following URL for progress report:
https://neos-server.org/neos/cgi-bin/nph-neos-solver.cgi?admin=results&jobnumber=17536350&pass=FaxWLjyN
Job 17536350 dispatched
password: FaxWLjyN
---------- Begin Solver Output -----------
Condor submit: 'neos.submit'
Condor submit: 'watchdog.submit'
Job submitted to NEOS HTCondor pool.
kestrel_options:solver=path
Supplied solver options are: threads=4
You are using the solver path.
Executing on prod-exec-6.neos-server.org
Path 4.7.03: threads=4
Path 4.7.03 (Sat May 31 16:42:27 2014)
Written by Todd Munson, Steven Dirkse, and Michael Ferris
INITIAL POINT STATISTICS
Maximum of X. . . . . . . . . . 0.0000e+00 var: (_svar[1])
Maximum of F. . . . . . . . . . 8.0000e+05 eqn: (_scon[6])
Maximum of Grad F . . . . . . . 1.3970e+03 eqn: (_scon[14])
var: (_svar[5])
INITIAL JACOBIAN NORM STATISTICS
Maximum Row Norm. . . . . . . . 7.8450e+03 eqn: (_scon[5])
Minimum Row Norm. . . . . . . . 1.2500e+02 eqn: (_scon[4])
Maximum Column Norm . . . . . . 7.8450e+03 var: (_svar[5])
Minimum Column Norm . . . . . . 1.2500e+02 var: (_svar[4])
Crash Log
major func diff size residual step prox (label)
0 0 1.8158e+06 0.0e+00 (_scon[6])
pn_search terminated: no progress.
Major Iteration Log
major minor func grad residual step type prox inorm (label)
0 0 13 1 1.8158e+06 I 1.0e+01 1.6e+06 (_scon[6])
1 17 14 2 1.1546e+04 1.0e+00 SO 4.0e+00 7.7e+03 (_scon[13])
2 7 15 3 1.1928e+03 1.0e+00 SO 1.6e+00 9.0e+02 (_scon[8])
3 5 16 4 3.7361e+02 1.0e+00 SO 6.4e-01 2.8e+02 (_scon[13])
4 1 17 5 2.6071e+00 1.0e+00 SO 2.6e-01 2.6e+00 (_scon[3])
5 1 18 6 1.3911e-02 1.0e+00 SO 1.0e-01 1.4e-02 (_scon[8])
6 1 19 7 1.1671e-06 1.0e+00 SO 1.4e-03 9.3e-07 (_scon[3])
7 1 20 8 1.4554e-11 1.0e+00 SO 1.2e-07 1.5e-11 (_scon[2])
FINAL STATISTICS
Inf-Norm of Complementarity . . 7.9248e-11 eqn: (_scon[2])
Inf-Norm of Normal Map. . . . . 1.4552e-11 eqn: (_scon[2])
Inf-Norm of Minimum Map . . . . 1.4552e-11 eqn: (_scon[2])
Inf-Norm of Fischer Function. . 1.4552e-11 eqn: (_scon[2])
Inf-Norm of Grad Fischer Fcn. . 5.8208e-10 eqn: (_scon[10])
Two-Norm of Grad Fischer Fcn. . 1.1748e-09
FINAL POINT STATISTICS
Maximum of X. . . . . . . . . . 1.8894e+03 var: (_svar[13])
Maximum of F. . . . . . . . . . 3.5116e+06 eqn: (_scon[5])
Maximum of Grad F . . . . . . . 1.3970e+03 eqn: (_scon[14])
var: (_svar[5])
** EXIT - solution found.
Major Iterations. . . . 7
Minor Iterations. . . . 33
Restarts. . . . . . . . 0
Crash Iterations. . . . 0
Gradient Steps. . . . . 0
Function Evaluations. . 20
Gradient Evaluations. . 8
Basis Time. . . . . . . 0.000151
Total Time. . . . . . . 0.000461
Residual. . . . . . . . 1.455369e-11
Path 4.7.03: Solution found.
7 iterations (0 for crash); 33 pivots.
20 function, 8 gradient evaluations.
Level [*] :=
P1 0
P1a 1555.3
P2 0
P2a 0
P2b 0
P3 147.465
P3c 1889.4
P4 0
;
sum{j in ACT} cost[j]*Level[j] = 6808640
assert ampl.solve_result == "solved", ampl.solve_result