{ "cells": [ { "cell_type": "markdown", "id": "3065afde", "metadata": {}, "source": [ "\n", "*This notebook contains material from [CBE60499](https://ndcbe.github.io/CBE60499);\n", "content is available [on Github](git@github.com:ndcbe/CBE60499.git).*\n" ] }, { "cell_type": "markdown", "id": "fa78f46f", "metadata": {}, "source": [ "\n", "< [4.5 Second Order Optimality Conditions](https://ndcbe.github.io/CBE60499/04.05-Second-Order.html) | [Contents](toc.html) | [Tag Index](tag_index.html) | [4.7 Simple Netwon Method for Equality Constrained NLPs](https://ndcbe.github.io/CBE60499/04.07-Interior-Point1.html) >
"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[4.6 NLP Diagnostics with Degeneracy Hunter](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6-NLP-Diagnostics-with-Degeneracy-Hunter)",
"section": "4.6 NLP Diagnostics with Degeneracy Hunter"
}
},
"source": [
"# 4.6 NLP Diagnostics with Degeneracy Hunter\n",
"\n",
"Created by Prof. Alex Dowling (adowling@nd.edu) at the University of Notre Dame. This notebook (and Degeneracy Hunter) are available through the [Institute for the Design of Advanced Energy Systems](https://idaes.org) under a BSD-3 license.\n",
"\n",
"This notebook shows how to use the following Degeneracy Hunter features using two motivating examples:\n",
"* Inspect constraint violations and bounds of a Pyomo model\n",
"* Compute the Irreducible Degenerate Set (IDS) for a Pyomo model\n",
"* Demonstrates the Ipopt performance benefits from removing a single redundant constraint\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbpages": {
"level": 1,
"link": "[4.6 NLP Diagnostics with Degeneracy Hunter](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6-NLP-Diagnostics-with-Degeneracy-Hunter)",
"section": "4.6 NLP Diagnostics with Degeneracy Hunter"
}
},
"outputs": [],
"source": [
"# This code cell installs packages on Colab\n",
"\n",
"import sys\n",
"if \"google.colab\" in sys.modules:\n",
" !wget \"https://raw.githubusercontent.com/ndcbe/CBE60499/main/notebooks/helper.py\"\n",
" import helper\n",
" helper.install_idaes()\n",
" helper.install_ipopt()\n",
" helper.install_cbc()"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.6.1 Setup](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.1-Setup)",
"section": "4.6.1 Setup"
}
},
"source": [
"## 4.6.1 Setup\n",
"\n",
"We start by importing Pyomo and Degeneracy Hunter."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.6.1 Setup](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.1-Setup)",
"section": "4.6.1 Setup"
}
},
"outputs": [],
"source": [
"import pyomo.environ as pyo\n",
"\n",
"from idaes.core.util.model_diagnostics import DegeneracyHunter\n",
"\n",
"milp_solver = pyo.SolverFactory('cbc')"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.6.2 Example 1: Well-Behaved Nonlinear Program](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2-Example-1:-Well-Behaved-Nonlinear-Program)",
"section": "4.6.2 Example 1: Well-Behaved Nonlinear Program"
}
},
"source": [
"## 4.6.2 Example 1: Well-Behaved Nonlinear Program"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.6.2 Example 1: Well-Behaved Nonlinear Program](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2-Example-1:-Well-Behaved-Nonlinear-Program)",
"section": "4.6.2 Example 1: Well-Behaved Nonlinear Program"
}
},
"source": [
"Consider the following \"well-behaved\" nonlinear optimization problem."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.6.2 Example 1: Well-Behaved Nonlinear Program](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2-Example-1:-Well-Behaved-Nonlinear-Program)",
"section": "4.6.2 Example 1: Well-Behaved Nonlinear Program"
}
},
"source": [
"$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{0,...,4\\}} x_i^2\\\\\n",
"\\mathrm{s.t.} \\quad & x_0 + x_1 - x_3 \\geq 10 \\\\\n",
"& x_0 \\times x_3 + x_1 \\geq 0 \\\\\n",
"& x_4 \\times x_3 + x_0 \\times x_3 + x_4 = 0\n",
"\\end{align*} $$"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[4.6.2 Example 1: Well-Behaved Nonlinear Program](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2-Example-1:-Well-Behaved-Nonlinear-Program)",
"section": "4.6.2 Example 1: Well-Behaved Nonlinear Program"
}
},
"source": [
"This problem is feasible, well-initialized, and standard constraint qualifications hold. As expected, we have no trouble solving this problem."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.1 Define the model in Pyomo](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.1-Define-the-model-in-Pyomo)",
"section": "4.6.2.1 Define the model in Pyomo"
}
},
"source": [
"### 4.6.2.1 Define the model in Pyomo"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.1 Define the model in Pyomo](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.1-Define-the-model-in-Pyomo)",
"section": "4.6.2.1 Define the model in Pyomo"
}
},
"source": [
"We start by defining the optimization problem in Pyomo."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.1 Define the model in Pyomo](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.1-Define-the-model-in-Pyomo)",
"section": "4.6.2.1 Define the model in Pyomo"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 Set Declarations\n",
" I : Size=1, Index=None, Ordered=Insertion\n",
" Key : Dimen : Domain : Size : Members\n",
" None : 1 : Any : 5 : {0, 1, 2, 3, 4}\n",
"\n",
"1 Var Declarations\n",
" x : Size=5, Index=I\n",
" Key : Lower : Value : Upper : Fixed : Stale : Domain\n",
" 0 : -10 : 1.0 : 10 : False : False : Reals\n",
" 1 : -10 : 1.0 : 10 : False : False : Reals\n",
" 2 : -10 : 1.0 : 10 : False : False : Reals\n",
" 3 : -10 : 1.0 : 10 : False : False : Reals\n",
" 4 : -10 : 1.0 : 10 : False : False : Reals\n",
"\n",
"1 Objective Declarations\n",
" obj : Size=1, Index=None, Active=True\n",
" Key : Active : Sense : Expression\n",
" None : True : minimize : x[0]**2 + x[1]**2 + x[2]**2 + x[3]**2 + x[4]**2\n",
"\n",
"3 Constraint Declarations\n",
" con1 : Size=1, Index=None, Active=True\n",
" Key : Lower : Body : Upper : Active\n",
" None : 10.0 : x[0] + x[1] - x[3] : +Inf : True\n",
" con2 : Size=1, Index=None, Active=True\n",
" Key : Lower : Body : Upper : Active\n",
" None : 0.0 : x[0]*x[3] + x[1] : +Inf : True\n",
" con3 : Size=1, Index=None, Active=True\n",
" Key : Lower : Body : Upper : Active\n",
" None : 0.0 : x[4]*x[3] + x[0]*x[3] - x[4] : 0.0 : True\n",
"\n",
"6 Declarations: I x con1 con2 con3 obj\n"
]
}
],
"source": [
"m = pyo.ConcreteModel()\n",
"\n",
"m.I = pyo.Set(initialize=[i for i in range(5)])\n",
"\n",
"m.x = pyo.Var(m.I,bounds=(-10,10),initialize=1.0)\n",
"\n",
"m.con1 = pyo.Constraint(expr=m.x[0] + m.x[1] - m.x[3] >= 10)\n",
"m.con2 = pyo.Constraint(expr=m.x[0]*m.x[3] + m.x[1] >= 0)\n",
"m.con3 = pyo.Constraint(expr=m.x[4]*m.x[3] + m.x[0]*m.x[3] - m.x[4] == 0)\n",
"\n",
"m.obj = pyo.Objective(expr=sum(m.x[i]**2 for i in m.I))\n",
"\n",
"m.pprint()"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.2 Evaluate the initial point](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.2-Evaluate-the-initial-point)",
"section": "4.6.2.2 Evaluate the initial point"
}
},
"source": [
"### 4.6.2.2 Evaluate the initial point"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.2 Evaluate the initial point](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.2-Evaluate-the-initial-point)",
"section": "4.6.2.2 Evaluate the initial point"
}
},
"source": [
"Initialization is extremely important for nonlinear optimization problems. By setting the Ipopt option `max_iter` to zero, we can inspect the initial point."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.2 Evaluate the initial point](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.2-Evaluate-the-initial-point)",
"section": "4.6.2.2 Evaluate the initial point"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ipopt 3.13.2: max_iter=0\n",
"\n",
"\n",
"******************************************************************************\n",
"This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
" Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
" For more information visit http://projects.coin-or.org/Ipopt\n",
"\n",
"This version of Ipopt was compiled from source code available at\n",
" https://github.com/IDAES/Ipopt as part of the Institute for the Design of\n",
" Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE\n",
" Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.\n",
"\n",
"This version of Ipopt was compiled using HSL, a collection of Fortran codes\n",
" for large-scale scientific computation. All technical papers, sales and\n",
" publicity material resulting from use of the HSL codes within IPOPT must\n",
" contain the following acknowledgement:\n",
" HSL, a collection of Fortran codes for large-scale scientific\n",
" computation. See http://www.hsl.rl.ac.uk.\n",
"******************************************************************************\n",
"\n",
"This is Ipopt version 3.13.2, running with linear solver ma27.\n",
"\n",
"Number of nonzeros in equality constraint Jacobian...: 3\n",
"Number of nonzeros in inequality constraint Jacobian.: 6\n",
"Number of nonzeros in Lagrangian Hessian.............: 7\n",
"\n",
"Total number of variables............................: 5\n",
" variables with only lower bounds: 0\n",
" variables with lower and upper bounds: 5\n",
" variables with only upper bounds: 0\n",
"Total number of equality constraints.................: 1\n",
"Total number of inequality constraints...............: 2\n",
" inequality constraints with only lower bounds: 2\n",
" inequality constraints with lower and upper bounds: 0\n",
" inequality constraints with only upper bounds: 0\n",
"\n",
"iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n",
" 0 5.0000000e+00 9.00e+00 2.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n",
"\n",
"Number of Iterations....: 0\n",
"\n",
" (scaled) (unscaled)\n",
"Objective...............: 5.0000000000000000e+00 5.0000000000000000e+00\n",
"Dual infeasibility......: 2.0000000000000000e+00 2.0000000000000000e+00\n",
"Constraint violation....: 8.9999999000000006e+00 8.9999999000000006e+00\n",
"Complementarity.........: 1.1000000099999999e+01 1.1000000099999999e+01\n",
"Overall NLP error.......: 1.1000000099999999e+01 1.1000000099999999e+01\n",
"\n",
"\n",
"Number of objective function evaluations = 1\n",
"Number of objective gradient evaluations = 1\n",
"Number of equality constraint evaluations = 1\n",
"Number of inequality constraint evaluations = 1\n",
"Number of equality constraint Jacobian evaluations = 1\n",
"Number of inequality constraint Jacobian evaluations = 1\n",
"Number of Lagrangian Hessian evaluations = 0\n",
"Total CPU secs in IPOPT (w/o function evaluations) = 0.000\n",
"Total CPU secs in NLP function evaluations = 0.000\n",
"\n",
"EXIT: Maximum Number of Iterations Exceeded.\n",
"WARNING: Loading a SolverResults object with a warning status into\n",
" model.name=\"unknown\";\n",
" - termination condition: maxIterations\n",
" - message from solver: Ipopt 3.13.2\\x3a Maximum Number of Iterations\n",
" Exceeded.\n"
]
}
],
"source": [
"# Specify Ipopt as the solver\n",
"opt = pyo.SolverFactory('ipopt')\n",
"\n",
"# Specifying an iteration limit of 0 allows us to inspect the initial point\n",
"opt.options['max_iter'] = 0\n",
"\n",
"# \"Solving\" the model with an iteration limit of 0 load the initial point and applies\n",
"# any preprocessors (e.g., enforces bounds)\n",
"opt.solve(m, tee=True)\n",
"\n",
"# Create Degeneracy Hunter object\n",
"dh = DegeneracyHunter(m, solver=milp_solver)"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.2 Evaluate the initial point](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.2-Evaluate-the-initial-point)",
"section": "4.6.2.2 Evaluate the initial point"
}
},
"source": [
"We expect the exit status `Maximum Number of Iterations Exceeded` because we told Ipopt to take zero iterations (only evaluate the initial point)."
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.3 Identify the constraint residuals larger than 0.1](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.3-Identify-the-constraint-residuals-larger-than-0.1)",
"section": "4.6.2.3 Identify the constraint residuals larger than 0.1"
}
},
"source": [
"### 4.6.2.3 Identify the constraint residuals larger than 0.1"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.3 Identify the constraint residuals larger than 0.1](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.3-Identify-the-constraint-residuals-larger-than-0.1)",
"section": "4.6.2.3 Identify the constraint residuals larger than 0.1"
}
},
"source": [
"When developing nonlinear optimization models, one often wants to know: \"what constraints are violated at the initial point (or more generally the point the solver terminated) within a given tolerance?\" Degeneracy Hunter makes this very easy by provided a simple interface to several IDAES utility functions.\n",
"\n",
"The following line of code will print out all constraints with residuals larger than `0.1`:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"nbpages": {
"level": 3,
"link": "[4.6.2.3 Identify the constraint residuals larger than 0.1](https://ndcbe.github.io/CBE60499/04.06-NLP-Diagnostics.html#4.6.2.3-Identify-the-constraint-residuals-larger-than-0.1)",
"section": "4.6.2.3 Identify the constraint residuals larger than 0.1"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" \n",
"All constraints with residuals larger than 0.1 :\n",
"\n",
"count = 0 \t|residual| = 9.0\n",
"con1 : Size=1, Index=None, Active=True\n",
" Key : Lower : Body : Upper : Active\n",
" None : 10.0 : x[0] + x[1] - x[3] : +Inf : True\n",
"variable\tlower\tvalue\tupper\n",
"x[0] \t\t -10 \t 1.0 \t 10\n",
"x[1] \t\t -10 \t 1.0 \t 10\n",
"x[3] \t\t -10 \t 1.0 \t 10\n",
"\n",
"count = 1 \t|residual| = 1.0\n",
"con3 : Size=1, Index=None, Active=True\n",
" Key : Lower : Body : Upper : Active\n",
" None : 0.0 : x[4]*x[3] + x[0]*x[3] - x[4] : 0.0 : True\n",
"variable\tlower\tvalue\tupper\n",
"x[4] \t\t -10 \t 1.0 \t 10\n",
"x[3] \t\t -10 \t 1.0 \t 10\n",
"x[0] \t\t -10 \t 1.0 \t 10\n",
"No constraints with residuals larger than 0.1 !\n"
]
},
{
"data": {
"text/plain": [
"dict_keys([
"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.8.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}