{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Simple nonlinear complementarity problem\n",
"\n",
"**Randall Romero Aguilar, PhD**\n",
"\n",
"This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.\n",
"\n",
"Original (Matlab) CompEcon file: **demslv13.m**\n",
"\n",
"Running this file requires the Python version of CompEcon. This can be installed with pip by running\n",
"\n",
" !pip install compecon --upgrade\n",
"\n",
"Last updated: 2022-Sept-05\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## About\n",
"\n",
"The problem is to solve\n",
"\n",
"\\begin{equation}\n",
"f(x, y) = \\begin{bmatrix}1+xy -2x^3-x\\\\ 2x^2-y\\end{bmatrix}\n",
"\\end{equation}\n",
"\n",
"subject to $0 \\leq x, y \\leq 1$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from compecon import MCP, jacobian\n",
"import numpy as np\n",
"\n",
"x0 = [0.5, 0.5] "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solving the problem without the Jacobian\n",
"To solve this problem we create a **MCP** object using a lambda function."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def func(z):\n",
" x, y = z\n",
" return [1 + x*y - 2*x**3 - x, 2*x**2 - y]\n",
" \n",
"F = MCP(func, [0, 0], [1,1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve for initial guess $x_0 = [0.5, 0.5]$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using the MINMAX transformation\n"
]
},
{
"ename": "AttributeError",
"evalue": "'list' object has no attribute 'size'",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mAttributeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_12340\\4017661448.py\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;31m#x = F.broyden(x0, transform='minmax', show=True)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mF\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbroyden\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtransform\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'ssmooth'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshow\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# FIXME: generates error\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 3\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34mf'Solution is {x=}.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\ProgramData\\Anaconda3\\lib\\site-packages\\compecon\\nonlinear.py\u001b[0m in \u001b[0;36mbroyden\u001b[1;34m(self, x0, **kwargs)\u001b[0m\n\u001b[0;32m 255\u001b[0m \u001b[0mmaxit\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmaxsteps\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mall_x\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mopts\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'maxit'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'maxsteps'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'tol'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'all_x'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 256\u001b[0m \u001b[0mmaxsteps\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmaxsteps\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 257\u001b[1;33m \u001b[0mfx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 258\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 259\u001b[0m \u001b[0mJinv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreset_inverse_jacobian\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\ProgramData\\Anaconda3\\lib\\site-packages\\compecon\\nonlinear.py\u001b[0m in \u001b[0;36m_minmax\u001b[1;34m(self, x)\u001b[0m\n\u001b[0;32m 507\u001b[0m \u001b[0mfx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mJ\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfx\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 508\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 509\u001b[1;33m \u001b[0mJ\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mjacobian\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_original\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 510\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 511\u001b[0m \u001b[0mfhat\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mda\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdb\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\ProgramData\\Anaconda3\\lib\\site-packages\\compecon\\tools.py\u001b[0m in \u001b[0;36mjacobian\u001b[1;34m(func, x, *args, **kwargs)\u001b[0m\n\u001b[0;32m 126\u001b[0m \u001b[0mdx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 127\u001b[0m \u001b[0mf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mF\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 128\u001b[1;33m \u001b[0mdf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 129\u001b[0m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 130\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mAttributeError\u001b[0m: 'list' object has no attribute 'size'"
]
}
],
"source": [
"#x = F.broyden(x0, transform='minmax', show=True) # FIXME: generates error\n",
"#x = F.broyden(x0, transform='ssmooth', show=True) # FIXME: generates error\n",
"#print(f'Solution is {x=}.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solving the problem with the Jacobian"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def func2(z):\n",
" x, y = z\n",
" f = [1 + x*y - 2*x**3 - x, 2*x**2 - y]\n",
" J = [[y-6*x**2-1, x],[4*x, -1]]\n",
" return f, J \n",
"\n",
"F2 = MCP(func2, [0, 0], [1,1])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"pycharm": {
"name": "#%%\n"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using the MINMAX transformation\n",
"Solving nonlinear equations by Newton's method\n",
"it bstep change\n",
"--------------------\n",
" 0 1 1.56e-01\n",
" 1 0 9.84e-03\n",
" 2 0 3.19e-05\n",
" 3 0 3.40e-10\n",
"Solution is x=array([0.7937, 1. ]).\n"
]
}
],
"source": [
"x = F2.zero(x0, transform='minmax', show=True)\n",
"print(f'Solution is {x=}.')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.9.7 ('base')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
},
"vscode": {
"interpreter": {
"hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
}
}
},
"nbformat": 4,
"nbformat_minor": 1
}
|