{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Nonlinear complementarity problem methods\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: **demslv08.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-04\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## About\n", "\n", "Solve nonlinear complementarity problem on $R^2$ using semismooth and minmax methods.\n", "\n", "Function to be solved is\n", "\n", "\\begin{equation*}\n", "f(x,y) = \\begin{bmatrix} 200x(y - x^2) + 1 - x\\\\100(x^2 - y)\\end{bmatrix}\n", "\\end{equation*}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set up the problem\n", "The class **MCP** is used to represent mixed-complementarity problems. To create one instance, we define the objective function and the boundaries $a$ and $b$ such that for $a \\leq x \\leq b$. Apart from the required parameters, we can specify options to be used when solving the problem." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from compecon import MCP, tic, toc" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def f(z):\n", " x, y = z\n", " fval = [200*x*(y - x**2) + 1 - x,\n", " 100*(x**2 - y)]\n", " fjac = [[200*(y - x**2) - 400*x**2 - 1, 200*x],\n", " [200*x, -100]]\n", "\n", " return fval, fjac" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generate problem test data" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "z = 2 * np.random.randn(2, 2)\n", "a = 1 + np.min(z, 0)\n", "b = 1 + np.max(z, 0)\n", "\n", "F = MCP(f, a, b, maxit=1500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solve by applying Newton method\n", "We'll use a random initial guess $x$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "F.x0 = np.random.randn(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Check the Jacobian " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "All numerical derivatives differ from\n", "the user-provided ones by less than 8 decimal digits.\n", "\n", "The maximum error is 7.30e-09, for row 0 and column 0.\n", "\n" ] } ], "source": [ "F.user_provides_jacobian = True\n", "F.check_jacobian()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using minmax formulation" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "t0 = tic()\n", "x1 = F.newton(transform='minmax')\n", "t1 = 100 * toc(t0)\n", "n1 = F.fnorm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using semismooth formulation" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "t0 = tic()\n", "x2 = F.newton(transform='ssmooth')\n", "t2 = 100*toc(t0)\n", "n2 = F.fnorm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Print results" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Hundreds of seconds required to solve nonlinear complementarity \n", "problem on R^2 using minmax and semismooth formulations, with \n", "randomly generated bounds \n", "\ta = [2.08, -0.08] \n", "\tb = [4.60, 4.42]\n" ] }, { "data": { "text/html": [ "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
 TimeNormx1x2
Algorithm    
Newton minmax1.1000.0000002.0834.337
Newton semismooth0.5000.0000002.0834.337
\n" ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(\n", " 'Hundreds of seconds required to solve nonlinear complementarity \\n' +\n", " 'problem on R^2 using minmax and semismooth formulations, with \\n' +\n", " 'randomly generated bounds \\n' +\n", " f'\\ta = [{a[0]:4.2f}, {a[1]:4.2f}] \\n' + \n", " f'\\tb = [{b[0]:4.2f}, {b[1]:4.2f}]'\n", " )\n", "\n", "results = pd.DataFrame.from_records(\n", " [('Newton minmax',t1, n1, *x1), \n", " ('Newton semismooth', t2, n2, *x2)],\n", " columns = ['Algorithm', 'Time', 'Norm', 'x1', 'x2']\n", " ).set_index('Algorithm')\n", "\n", "results.style.format(precision=3, subset=['Time', 'x1', 'x2'])" ] } ], "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": 0 }