diff --git a/.github/actions/install/action.yml b/.github/actions/install/action.yml index c3100786..68e497cf 100644 --- a/.github/actions/install/action.yml +++ b/.github/actions/install/action.yml @@ -42,3 +42,9 @@ runs: echo '::group::Output of "idaes get-extensions" command' idaes get-extensions --extra petsc --verbose echo '::endgroup::' + - name: Install SCIP from AMPL + shell: bash -l {0} + run: | + echo '::group::Output of "pip install ampl_module_scip" command' + ${{ inputs.install-command }} --index-url https://pypi.ampl.com ampl_module_scip + echo '::endgroup::' \ No newline at end of file diff --git a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter.ipynb b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter.ipynb index 063eb4b2..d09a2a4a 100644 --- a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter.ipynb +++ b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": { "tags": [ "header", @@ -34,324 +34,97 @@ "Maintainer: Andrew Lee \n", "Updated: 2023-06-01 \n", "\n", - "This notebook shows how to use the following Degeneracy Hunter features using two motivating examples:\n", + "> Alexander W. Dowling and Lorenz T. Biegler (2015). Degeneracy Hunter: An Algorithm for Determining Irreducible Sets of Degenerate Constraints in Mathematical Programs. 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering. Ed. by Krist V. Gernaey, Jakob K. Huusom, and Raqul Gani. Computer-Aided Chemical Engineering, Vol. 37, p. 809 – 814. [link to PDF](https://dowlinglab.nd.edu/assets/447116/2015_degeneracy_hunter_an_algorithm_for_determining_irreducible_sets_of_degenerate_constraints_in_mathematical_prog.pdf)\n", + "\n", + "This notebook shows how to use the Degeneracy Hunter tool which is part of the IDAES Diagnostics Toolbox features using two motivating examples:\n", + "\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": "markdown", - "metadata": {}, - "source": [ - "## Setup\n", + "* Demonstrates the performance benefits in Ipopt from removing a single redundant constraint\n", "\n", - "We start by importing Pyomo and Degeneracy Hunter." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import pyomo.environ as pyo\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver. This example uses the open-source solver SCIP (https://scipopt.org/), however users may specify another solver of their choice if they have one available.\n", + "
\n", "\n", - "from idaes.core.util.model_diagnostics import DegeneracyHunter" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example 1: Well-Behaved Nonlinear Program" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Consider the following \"well-behaved\" nonlinear optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "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": {}, - "source": [ - "### Define the model in Pyomo" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We start by defining the optimization problem in Pyomo." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "m = pyo.ConcreteModel()\n", + "## What does Degeneracy Mean?\n", "\n", - "m.I = pyo.Set(initialize=[i for i in range(5)])\n", + "In simple terms, a degenerate model contains one or more constraints that do not add any additional information that is not already included in other constraints in the model. The simplest example of this is the case where a constraint in the model is duplicated. For example, consider the model below:\n", "\n", - "m.x = pyo.Var(m.I, bounds=(-10, 10), initialize=1.0)\n", + "$$\\begin{align*}\n", + "& x_1 + x_2 = 1 \\\\\n", + "& x_1 + x_2 = 1 \\\\\n", + "\\end{align*} $$\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", + "Notice the two equality constraints are identical and thus redundant. This means the constraint qualifications that solvers rely on (e.g., [Linear Independence Constraint Qualification or LICQ](https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions#Regularity_conditions_(or_constraint_qualifications))) do not hold which has three important implications:\n", "\n", - "m.obj = pyo.Objective(expr=sum(m.x[i] ** 2 for i in m.I))\n", + "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", + "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", + "3. Theoretical convergence properties of optimization algorithms may not hold\n", "\n", - "m.pprint()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Evaluate the initial point" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": 3, - "metadata": {}, - "outputs": [], - "source": [ - "# Specify Ipopt as the solver\n", - "opt = pyo.SolverFactory(\"ipopt\")\n", + "Put another way, for a deterministic (square) problem we required that the number of equality constraints be equal to the number of free variables. In this case, there appear to be two equality constraints but in reality there is effectively only one as the two constraints are the same. Thus, the problem is not well defined and there exist a family of solutions that can satisfy the constraints given; any combination of values for ``x1`` and ``x2`` that sum to one will satisfy the constraints as written.\n", "\n", - "# Specifying an iteration limit of 0 allows us to inspect the initial point\n", - "opt.options[\"max_iter\"] = 0\n", + "In practice, degeneracies are rarely as straight-forward as a duplicated constraint. A more common occurrence is a set of linearly dependent constraints; that is a system of constraints where one constraint in the set can be formed by a linear combination of other constraints in the set. For example:\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", + "$$\\begin{align*}\n", + "& x_1 + x_2 + x_3 = 1 \\\\\n", + "& 2x_1 + 2x_2 + 2x_3 = 2 \\\\\n", + "\\end{align*} $$\n", "\n", - "# Create Degeneracy Hunter object\n", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh = DegeneracyHunter(m, solver=pyo.SolverFactory(\"cbc\"))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "source": [ - "### Identify the constraint residuals larger than 0.1" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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", + "Here, the second constraint can be formed by multiplying the first constraint by two. Another common example of linearly dependent constraints is the combination of a set of material balances for all components in a system and an overall material balance. In this case, the overall material balances can be formed by adding all the individual component balances together. Depending on the number of components being modeled, this system of constraints can be quite large, but the end result is the same; the model effectively has one less constraint than it would appear to have.\n", "\n", - "The following line of code will print out all constraints with residuals larger than `0.1`:" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Important: Ipopt does several preprocessing steps when we executed it with zero iterations. When checking the initial point, it is strongly recommended to call Ipopt with zero iterations first. Otherwise, you will not be analyzing the initial point Ipopt starts with." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify all variables within 1 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Another common question when developing optimization models is, \"Which variables are within their bounds by a given tolerance?\" Below is the syntax:" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_variable_bounds(tol=1.0)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve the optimization problem" + "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the purpose of Degeneracy Hunter. Let's see it in action." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now we can solve the optimization problem. We first set the number of iterations to 50 and then resolve with Ipopt." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m, tee=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, Ipopt has no trouble solving this optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check if any constraint residuals are large than 1E-14" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's now inspect the new solution to see which (if any) constraints have residuals larger than 10$^{-14}$." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_residuals(tol=1e-14)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, all of the constraints are satisfied, even with this fairly tight tolerance." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify all variables within 1E-5 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, let's check if any of the variables are near their bounds at the new solution." + "## Setup\n", + "\n", + "We start by importing Pyomo and the IDAES Diagnostics Toolbox." ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "dh.check_variable_bounds(tol=1e-5)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Great, no variables are near their bounds. If a variable was at its bound, it is important the inspect the model and confirm the bound is physically sensible/what you intended." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check the rank of the constraint Jacobian at the solution" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The main feature of Degeneracy Hunter is to check if an optimization problem is poorly formulated. Let's see what happens when we check the rank on a carefully formulated optimization problem:" + "import pyomo.environ as pyo\n", + "\n", + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox" ] }, { "cell_type": "code", - "execution_count": 9, - "metadata": {}, + "execution_count": 3, + "metadata": { + "tags": [ + "testing" + ] + }, "outputs": [], "source": [ - "dh.check_rank_equality_constraints()" + "from idaes.core.util.testing import _enable_scip_solver_for_testing\n", + "scip_solver = pyo.SolverFactory(\"scip\")\n", + "if not scip_solver.available():\n", + " _enable_scip_solver_for_testing()\n", + "assert scip_solver.available()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Example 2: Linear Program with Redundant Equality Constraints" + "## Example: Linear Program with Redundant Equality Constraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's apply Degeneracy Hunter to a poorly formulated optimization problem:\n", + "Let's apply these tools to a poorly formulated optimization problem:\n", "\n", "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", @@ -361,256 +134,661 @@ "& x_1 + x_2 + x_3 = 1 \\\\\n", "\\end{align*} $$\n", "\n", + "As you can see, this is similar to our example above with duplicated equality constraints.\n", "\n", - "Notice the two equality constraints are redundant. This means the constraint qualifications (e.g., LICQ) do not hold which has three important implications:\n", - "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", - "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", - "3. Theoretical convergence properties of optimization algorithms may not hold\n", - "\n", - "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the primary purpose of Degeneracy Hunter. Let's see it in action." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Define the model in Pyomo" + "The cell below shows how to construct this model in Pyomo." ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ - "def example2(with_degenerate_constraint=True):\n", - " \"\"\"Create the Pyomo model for Example 2\n", - "\n", - " Arguments:\n", - " with_degenerate_constraint: Boolean, if True, include the redundant linear constraint\n", - "\n", - " Returns:\n", - " m2: Pyomo model\n", - " \"\"\"\n", - "\n", - " m2 = pyo.ConcreteModel()\n", - "\n", - " m2.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", - "\n", - " m2.x = pyo.Var(m2.I, bounds=(0, 5), initialize=1.0)\n", - "\n", - " m2.con1 = pyo.Constraint(expr=m2.x[1] + m2.x[2] >= 1)\n", - " m2.con2 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", - " m2.con3 = pyo.Constraint(expr=m2.x[2] - 2 * m2.x[3] <= 1)\n", - " m2.con4 = pyo.Constraint(expr=m2.x[1] + m2.x[3] >= 1)\n", - "\n", - " if with_degenerate_constraint:\n", - " m2.con5 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", + "def build_example(include_redundant_constraint=True):\n", + " m = pyo.ConcreteModel()\n", "\n", - " m2.obj = pyo.Objective(expr=sum(m2.x[i] for i in m2.I))\n", + " m.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", "\n", - " m2.pprint()\n", + " m.x = pyo.Var(m.I, bounds=(0, 5), initialize=1.0)\n", "\n", - " return m2\n", + " m.con1 = pyo.Constraint(expr=m.x[1] + m.x[2] >= 1)\n", + " m.con2 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", + " m.con3 = pyo.Constraint(expr=m.x[2] - 2 * m.x[3] <= 1)\n", + " m.con4 = pyo.Constraint(expr=m.x[1] + m.x[3] >= 1)\n", + " \n", + " if include_redundant_constraint:\n", + " m.con5 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", "\n", + " m.obj = pyo.Objective(expr=sum(m.x[i] for i in m.I))\n", + " \n", + " return m\n", "\n", - "# Create the Pyomo model for Example 2 including the redundant constraint\n", - "m2 = example2()" + "m = build_example()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Evaluate the initial point" + "### Check for Structural Issues\n", + "\n", + "Before we try to solve the model, the first thing we should do is use the Diagnostics Toolbox to check for any structural issues in our model. The cell below creates an instance of the Diagnostics Toolbox and calls the ``report_structural_issues`` method." ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 5, "metadata": {}, - "outputs": [], - "source": [ + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "Model Statistics\n", + "\n", + " Activated Blocks: 1 (Deactivated: 0)\n", + " Free Variables in Activated Constraints: 3 (External: 0)\n", + " Free Variables with only lower bounds: 0\n", + " Free Variables with only upper bounds: 0\n", + " Free Variables with upper and lower bounds: 3\n", + " Fixed Variables in Activated Constraints: 0 (External: 0)\n", + " Activated Equality Constraints: 2 (Deactivated: 0)\n", + " Activated Inequality Constraints: 3 (Deactivated: 0)\n", + " Activated Objectives: 1 (Deactivated: 0)\n", + "\n", + "------------------------------------------------------------------------------------\n", + "2 WARNINGS\n", + "\n", + " WARNING: 1 Degree of Freedom\n", + " WARNING: Structural singularity found\n", + " Under-Constrained Set: 3 variables, 2 constraints\n", + " Over-Constrained Set: 0 variables, 0 constraints\n", + "\n", + "------------------------------------------------------------------------------------\n", + "0 Cautions\n", + "\n", + " No cautions found!\n", + "\n", + "------------------------------------------------------------------------------------\n", + "Suggested next steps:\n", + "\n", + " display_underconstrained_set()\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "dt = DiagnosticsToolbox(m)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, the toolbox is reporting one degree of freedom and a structural singularity. In this case, these are both expected as we are trying to solve an optimization problem with one degree of freedom so we can move on to the next step. \n", + "\n", + "### Evaluate the initial point\n", + "\n", + "Before we try to solve our degenerate model, it can often be useful to check the state of the model at the initial step seen by the solver. The cell below creates an instance of Ipopt and sets the iteration limit to 0 so that it terminates at the initial condition." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.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...............: 3.0000000000000000e+00 3.0000000000000000e+00\n", + "Dual infeasibility......: 1.0000000000000000e+00 1.0000000000000000e+00\n", + "Constraint violation....: 2.0000000000000000e+00 2.0000000000000000e+00\n", + "Complementarity.........: 4.0000000499999997e+00 4.0000000499999997e+00\n", + "Overall NLP error.......: 4.0000000499999997e+00 4.0000000499999997e+00\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" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'warning', 'Message': 'Ipopt 3.13.2\\\\x3a Maximum Number of Iterations Exceeded.', 'Termination condition': 'maxIterations', 'Id': 400, 'Error rc': 0, 'Time': 0.004778385162353516}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver = 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", + "solver.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(m2, tee=True)\n", - "\n", - "# Create Degeneracy Hunter object\n", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh2 = DegeneracyHunter(m2, solver=pyo.SolverFactory(\"cbc\"))" + "solver.solve(m, tee=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Identify constraints with residuals greater than 0.1 at the initial point" + "### Identify constraints with large residuals at the initial point\n", + "\n", + "With the solution at the initial point, we can then use the Diagnostics Toolbox to check the residuals of all constraints at the initial point." ] }, { "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "dh2.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 7, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following constraint(s) have large residuals (>1.0E-05):\n", + "\n", + " con2: 2.00000E+00\n", + " con5: 2.00000E+00\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "### Solve the optimization problem and extract the solution" + "# Check for large residuals\n", + "dt.display_constraints_with_large_residuals()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's solve the optimization problem." + "We can also check for variables which are close to (or outside of) their bounds." ] }, { "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2, tee=True)\n", - "\n", - "for i in m2.I:\n", - " print(\"x[\", i, \"]=\", m2.x[i]())" - ] - }, - { - "cell_type": "markdown", + "execution_count": 8, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following variable(s) have values close to their bounds (abs=1.0E-04, rel=1.0E-04):\n", + "\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of linesearch evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution." + "dt.display_variables_near_bounds()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Check the rank of the Jacobian of the equality constraints" + "### Solve the optimization problem and extract the solution\n", + "\n", + "There appear to be no significant issues at the initial point, so let us move on to solving the optimization problem." ] }, { "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "n_deficient = dh2.check_rank_equality_constraints()" - ] - }, - { - "cell_type": "markdown", + "execution_count": 9, "metadata": {}, - "source": [ - "A singular value near 0 indicates the Jacobian of the equality constraints is rank deficient. For each near-zero singular value, there is likely one degenerate constraint." + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 2.17e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01h 1\n", + " 2 1.0184419e+00 1.84e-02 1.26e-02 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.82e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 4.99e-04 - 1.00e+00 5.82e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.33e+01 -2.5 1.88e-04 - 1.00e+00 2.94e-01f 2\n", + " 8 1.0000153e+00 1.53e-05 2.19e+01 -2.5 6.85e-05 - 1.00e+00 8.08e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.78e+02 -2.5 3.47e-05 - 1.00e+00 2.45e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000014e+00 1.44e-06 2.83e-08 -2.5 7.79e-06 - 1.00e+00 1.00e+00h 1\n", + " 11 1.0000000e+00 1.39e-09 1.84e-11 -5.7 2.56e-06 - 1.00e+00 1.00e+00h 1\n", + " 12 1.0000000e+00 1.39e-09 1.24e+02 -8.6 1.16e-08 - 1.00e+00 9.54e-07h 21\n", + " 13 9.9999999e-01 9.82e-09 1.14e-13 -8.6 1.33e-08 - 1.00e+00 1.00e+00s 22\n", + " 14 9.9999999e-01 9.96e-09 6.46e-14 -8.6 2.29e-10 - 1.00e+00 1.00e+00s 22\n", + " 15 9.9999999e-01 9.96e-09 1.82e+02 -9.0 4.00e-11 - 1.00e+00 9.54e-07h 21\n", + " 16 9.9999999e-01 1.00e-08 1.52e-13 -9.0 5.89e-11 - 1.00e+00 1.00e+00s 22\n", + "\n", + "Number of Iterations....: 16\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999000238915e-01 9.9999999000238915e-01\n", + "Dual infeasibility......: 1.5188693260016487e-13 1.5188693260016487e-13\n", + "Constraint violation....: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "Complementarity.........: 9.2217032157601989e-10 9.2217032157601989e-10\n", + "Overall NLP error.......: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 111\n", + "Number of objective gradient evaluations = 17\n", + "Number of equality constraint evaluations = 111\n", + "Number of inequality constraint evaluations = 111\n", + "Number of equality constraint Jacobian evaluations = 17\n", + "Number of inequality constraint Jacobian evaluations = 17\n", + "Number of Lagrangian Hessian evaluations = 16\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: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.04887509346008301}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver.options[\"max_iter\"] = 50\n", + "solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of line search evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution.\n", + "\n", + "The cell below shows the values of the variables at the solution." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 10, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000099975996 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 0.0 : 5 : False : False : Reals\n" + ] + } + ], "source": [ - "### Identify candidate degenerate constraints" + "m.x.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Degeneracy Hunter first identifies candidate degenerate equations." + "### Check the rank of the Jacobian of the equality constraints\n", + "\n", + "One way to check for the presence of degenerate equations is to calculate the rank of the Jacobian matrix of the equality constraints. A singular value near 0 indicates the Jacobian is rank deficient, and for each near-zero singular value there is likely one degenerate constraint.\n", + "\n", + "The Diagnostics Toolbox contains a Singular Value Decomposition (SVD) toolbox which can be used to determine the number of small singular values in the Jacobian for a model as shown below." ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 11, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "\n", + "Number of Singular Values less than 1.0E-06 is 1\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "ds2 = dh2.find_candidate_equations(verbose=True, tee=True)" + "svd = dt.prepare_svd_toolbox()\n", + "svd.display_rank_of_equality_constraints()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Find irreducible degenerate sets (IDS)" + "As can be seen above, there is one singular value less than 1e-6, suggesting we have one degenerate constraint." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Next, Degeneracy Hunter enumerates through the candidate equations. For each candidate equation, Degenerate Hunter solves a MILP to compute the irreducible degenerate set that must contain the candidate equation." + "### Finding Irreducible Degenerate Sets (IDS)\n", + "\n", + "Next, we can use Degeneracy Hunter to identify irreducible degenerate sets; that is, the smallest set of constraints that contain a redundant constraint.\n", + "\n", + "Degeneracy Hunter first identifies candidate degenerate equations by solving a mixed integer linear program (MILP). It then iterates through the candidate equations and solves a second MILP for each to compute the irreducible degenerate set that must contain the candidate equation.\n", + "\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver, users can specify their solver of choice via a configuration argument. The default option is SCIP (https://scipopt.org/) which is an open-source solver.\n", + "
\n", + "\n", + "The cell below shows how to create an instance of Degeneracy Hunter and generate a report of the irreducible degenerate sets." ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 12, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Candidate Equations\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving Candidates MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Irreducible Degenerate Sets\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model to compute irreducible degenerate set.\n", + "Solving MILP 1 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con2.\n", + "Solving MILP 2 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con5.\n", + "====================================================================================\n", + "Irreducible Degenerate Sets\n", + "\n", + " Irreducible Degenerate Set 0\n", + " nu Constraint Name\n", + " 1.0 con2\n", + " -1.0 con5\n", + "\n", + " Irreducible Degenerate Set 1\n", + " nu Constraint Name\n", + " -1.0 con2\n", + " 1.0 con5\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "ids = dh2.find_irreducible_degenerate_sets(verbose=True)" + "dh = dt.prepare_degeneracy_hunter()\n", + "dh.report_irreducible_degenerate_sets()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Reformulate Example 2" + "As can be seen above, Degeneracy Hunter identified two constraints as potentially redundant (the duplicate constraints ``con2`` and ``con5``). Degeneracy Hunter then reports an Irreducible Degenerate Set for each of these constraints, which in the case as identical. More complex models may have partially overlapping IDS's.\n", + "\n", + "These results show us that ``con2`` and ``con5`` are mutually redundant; i.e., we can eliminate one of these with no effect on the model solution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's reformulate the model by skipping/removing the redundant equality constraint:\n", + "### Reformulate Model to Remove Redundant Constraint\n", "\n", - "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", - "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", - "& x_1 + x_2 + x_3 = 1 \\\\\n", - "& x_2 - 2 x_3 \\leq 1 \\\\\n", - "& x_1 + x_3 \\geq 1\n", - "\\end{align*} $$" + "Now let's reformulate the model by removing the redundant equality constraint ``con5`` and solve the reformulated model." ] }, { "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [], - "source": [ - "m2b = example2(with_degenerate_constraint=False)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 13, "metadata": {}, - "source": [ - "### Solve the reformulated model" + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 1\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 6.30e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 1.34e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01f 1\n", + " 2 1.0184419e+00 1.84e-02 9.15e-03 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.88e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 5.00e-04 - 1.00e+00 5.81e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.34e+01 -2.5 1.88e-04 - 1.00e+00 2.93e-01f 2\n", + " 8 1.0000154e+00 1.54e-05 2.26e+01 -2.5 6.89e-05 - 1.00e+00 8.04e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.98e+02 -2.5 3.64e-05 - 1.00e+00 2.33e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000020e+00 2.04e-06 1.33e+02 -2.5 9.72e-06 - 1.00e+00 8.28e-01h 1\n", + " 11 1.0000016e+00 1.61e-06 2.26e+03 -2.5 4.79e-06 - 1.00e+00 2.12e-01f 2\n", + " 12 1.0000002e+00 2.46e-07 8.72e+02 -2.5 1.20e-06 - 1.00e+00 8.47e-01h 1\n", + " 13 1.0000002e+00 1.95e-07 1.71e+04 -2.5 7.02e-07 - 1.00e+00 2.09e-01f 2\n", + " 14 1.0000000e+00 1.54e-08 3.23e+03 -2.5 1.50e-07 - 1.00e+00 9.21e-01h 1\n", + " 15 1.0000000e+00 1.15e-08 9.99e+04 -2.5 6.89e-08 - 1.00e+00 2.54e-01f 2\n", + " 16 1.0000000e+00 2.22e-16 2.83e-08 -2.5 8.21e-09 - 1.00e+00 1.00e+00h 1\n", + " 17 1.0000000e+00 2.22e-16 4.14e-11 -8.6 8.25e-16 - 1.00e+00 1.00e+00 0\n", + "\n", + "Number of Iterations....: 17\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999999999978e-01 9.9999999999999978e-01\n", + "Dual infeasibility......: 4.1425156707686206e-11 4.1425156707686206e-11\n", + "Constraint violation....: 2.2204460492503131e-16 2.2204460492503131e-16\n", + "Complementarity.........: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "Overall NLP error.......: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 23\n", + "Number of objective gradient evaluations = 18\n", + "Number of equality constraint evaluations = 23\n", + "Number of inequality constraint evaluations = 23\n", + "Number of equality constraint Jacobian evaluations = 18\n", + "Number of inequality constraint Jacobian evaluations = 18\n", + "Number of Lagrangian Hessian evaluations = 17\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.002\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 4, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.006506204605102539}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m2 = build_example(include_redundant_constraint=False)\n", + "\n", + "solver.solve(m2, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let us check to see if the optimal solution has changed." ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000005924994 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 4.55844318120227e-11 : 5 : False : False : Reals\n" + ] + } + ], "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2b, tee=True)\n", - "\n", - "for i in m2b.I:\n", - " print(\"x[\", i, \"]=\", m.x[i]())" + "m2.x.display()" ] }, { @@ -650,13 +828,23 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, + "execution_count": 15, + "metadata": { + "tags": [ + "testing" + ] + }, "outputs": [], - "source": [] + "source": [ + "import pytest\n", + "assert pyo.value(m2.x[1]) == pytest.approx(1, rel=1e-6)\n", + "assert pyo.value(m2.x[2]) == pytest.approx(0, abs=1e-6)\n", + "assert pyo.value(m2.x[3]) == pytest.approx(0, abs=1e-6)" + ] } ], "metadata": { + "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", @@ -672,7 +860,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.6" + "version": "3.11.5" } }, "nbformat": 4, diff --git a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_doc.ipynb b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_doc.ipynb index de97c522..ac2388b6 100644 --- a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_doc.ipynb +++ b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_doc.ipynb @@ -1,1638 +1,835 @@ { - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": { - "tags": [ - "header", - "hide-cell" - ] - }, - "outputs": [], - "source": [ - "###############################################################################\n", - "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", - "# Framework (IDAES IP) was produced under the DOE Institute for the\n", - "# Design of Advanced Energy Systems (IDAES).\n", - "#\n", - "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", - "# University of California, through Lawrence Berkeley National Laboratory,\n", - "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", - "# University, West Virginia University Research Corporation, et al.\n", - "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", - "# for full copyright and license information.\n", - "###############################################################################" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Degeneracy Hunter Examples\n", - "Author: Prof. Alex Dowling (adowling@nd.edu), University of Notre Dame \n", - "Maintainer: Andrew Lee \n", - "Updated: 2023-06-01 \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": "markdown", - "metadata": {}, - "source": [ - "## Setup\n", - "\n", - "We start by importing Pyomo and Degeneracy Hunter." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "import pyomo.environ as pyo\n", - "\n", - "from idaes.core.util.model_diagnostics import DegeneracyHunter" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example 1: Well-Behaved Nonlinear Program" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Consider the following \"well-behaved\" nonlinear optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "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": {}, - "source": [ - "### Define the model in Pyomo" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We start by defining the optimization problem in Pyomo." - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "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": {}, - "source": [ - "### Evaluate the initial point" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "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" - ] + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "tags": [ + "header", + "hide-cell" + ] + }, + "outputs": [], + "source": [ + "###############################################################################\n", + "# The Institute for the Design of Advanced Energy Systems Integrated Platform\n", + "# Framework (IDAES IP) was produced under the DOE Institute for the\n", + "# Design of Advanced Energy Systems (IDAES).\n", + "#\n", + "# Copyright (c) 2018-2023 by the software owners: The Regents of the\n", + "# University of California, through Lawrence Berkeley National Laboratory,\n", + "# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon\n", + "# University, West Virginia University Research Corporation, et al.\n", + "# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md\n", + "# for full copyright and license information.\n", + "###############################################################################" + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "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" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Degeneracy Hunter Examples\n", + "Author: Prof. Alex Dowling (adowling@nd.edu), University of Notre Dame \n", + "Maintainer: Andrew Lee \n", + "Updated: 2023-06-01 \n", + "\n", + "> Alexander W. Dowling and Lorenz T. Biegler (2015). Degeneracy Hunter: An Algorithm for Determining Irreducible Sets of Degenerate Constraints in Mathematical Programs. 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering. Ed. by Krist V. Gernaey, Jakob K. Huusom, and Raqul Gani. Computer-Aided Chemical Engineering, Vol. 37, p. 809 \u2013 814. [link to PDF](https://dowlinglab.nd.edu/assets/447116/2015_degeneracy_hunter_an_algorithm_for_determining_irreducible_sets_of_degenerate_constraints_in_mathematical_prog.pdf)\n", + "\n", + "This notebook shows how to use the Degeneracy Hunter tool which is part of the IDAES Diagnostics Toolbox features using two motivating examples:\n", + "\n", + "* Inspect constraint violations and bounds of a Pyomo model\n", + "* Compute the Irreducible Degenerate Set (IDS) for a Pyomo model\n", + "* Demonstrates the performance benefits in Ipopt from removing a single redundant constraint\n", + "\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver. This example uses the open-source solver SCIP (https://scipopt.org/), however users may specify another solver of their choice if they have one available.\n", + "
\n", + "\n", + "## What does Degeneracy Mean?\n", + "\n", + "In simple terms, a degenerate model contains one or more constraints that do not add any additional information that is not already included in other constraints in the model. The simplest example of this is the case where a constraint in the model is duplicated. For example, consider the model below:\n", + "\n", + "$$\\begin{align*}\n", + "& x_1 + x_2 = 1 \\\\\n", + "& x_1 + x_2 = 1 \\\\\n", + "\\end{align*} $$\n", + "\n", + "Notice the two equality constraints are identical and thus redundant. This means the constraint qualifications that solvers rely on (e.g., [Linear Independence Constraint Qualification or LICQ](https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions#Regularity_conditions_(or_constraint_qualifications))) do not hold which has three important implications:\n", + "\n", + "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", + "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", + "3. Theoretical convergence properties of optimization algorithms may not hold\n", + "\n", + "Put another way, for a deterministic (square) problem we required that the number of equality constraints be equal to the number of free variables. In this case, there appear to be two equality constraints but in reality there is effectively only one as the two constraints are the same. Thus, the problem is not well defined and there exist a family of solutions that can satisfy the constraints given; any combination of values for ``x1`` and ``x2`` that sum to one will satisfy the constraints as written.\n", + "\n", + "In practice, degeneracies are rarely as straight-forward as a duplicated constraint. A more common occurrence is a set of linearly dependent constraints; that is a system of constraints where one constraint in the set can be formed by a linear combination of other constraints in the set. For example:\n", + "\n", + "$$\\begin{align*}\n", + "& x_1 + x_2 + x_3 = 1 \\\\\n", + "& 2x_1 + 2x_2 + 2x_3 = 2 \\\\\n", + "\\end{align*} $$\n", + "\n", + "Here, the second constraint can be formed by multiplying the first constraint by two. Another common example of linearly dependent constraints is the combination of a set of material balances for all components in a system and an overall material balance. In this case, the overall material balances can be formed by adding all the individual component balances together. Depending on the number of components being modeled, this system of constraints can be quite large, but the end result is the same; the model effectively has one less constraint than it would appear to have.\n", + "\n", + "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the purpose of Degeneracy Hunter. Let's see it in action." + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "2023-11-02 10:25:10 [WARNING] idaes.core.util.model_diagnostics: DEPRECATED: DegeneracyHunter is being deprecated in favor of the new\n", - "DiagnosticsToolbox. (deprecated in 2.2.0, will be removed in (or\n", - "after) 3.0.0)\n", - "(called from C:\\Users\\dkgun\\AppData\\Local\\Temp\\ipykernel_30328\\560750228.py:14)\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", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh = DegeneracyHunter(m, solver=pyo.SolverFactory(\"cbc\"))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "source": [ - "### Identify the constraint residuals larger than 0.1" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " \n", - "All constraints with residuals larger than 0.1 :\n", - "Count\tName\t|residual|\n", - "0 \t con1 \t 9.0\n", - "1 \t con3 \t 1.0\n" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "We start by importing Pyomo and the IDAES Diagnostics Toolbox." + ] }, { - "data": { - "text/plain": [ - "dict_keys([, ])" + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import pyomo.environ as pyo\n", + "\n", + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox" ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dh.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Important: Ipopt does several preprocessing steps when we executed it with zero iterations. When checking the initial point, it is strongly recommended to call Ipopt with zero iterations first. Otherwise, you will not be analyzing the initial point Ipopt starts with." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify all variables within 1 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Another common question when developing optimization models is, \"Which variables are within their bounds by a given tolerance?\" Below is the syntax:" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "2023-11-02 10:25:10 [WARNING] idaes.core.util.model_statistics: DEPRECATED: variables_near_bounds_generator has deprecated the\n", - "relative argument. Please set abs_tol and rel_tol arguments instead.\n", - "(deprecated in 2.2.0, will be removed in (or after) 3.0.0)\n", - "(called from C:\\Users\\dkgun\\miniconda3\\envs\\idaes-examples-py311\\Lib\\site-packages\\pyomo\\common\\collections\\component_set.py:61)\n" - ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "2023-11-02 10:25:10 [WARNING] idaes.core.util.model_statistics: DEPRECATED: variables_near_bounds_generator has deprecated the tol\n", - "argument. Please set abs_tol and rel_tol arguments instead.\n", - "(deprecated in 2.2.0, will be removed in (or after) 3.0.0)\n", - "(called from C:\\Users\\dkgun\\miniconda3\\envs\\idaes-examples-py311\\Lib\\site-packages\\pyomo\\common\\collections\\component_set.py:61)\n" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example: Linear Program with Redundant Equality Constraints" + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - " \n", - "Variables within 1.0 (absolute) of their bounds:\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[2] \t\t -10 \t 1.0 \t 10\n", - "x[3] \t\t -10 \t 1.0 \t 10\n", - "x[4] \t\t -10 \t 1.0 \t 10\n" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's apply these tools to a poorly formulated optimization problem:\n", + "\n", + "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", + "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", + "& x_1 + x_2 + x_3 = 1 \\\\\n", + "& x_2 - 2 x_3 \\leq 1 \\\\\n", + "& x_1 + x_3 \\geq 1 \\\\\n", + "& x_1 + x_2 + x_3 = 1 \\\\\n", + "\\end{align*} $$\n", + "\n", + "As you can see, this is similar to our example above with duplicated equality constraints.\n", + "\n", + "The cell below shows how to construct this model in Pyomo." + ] }, { - "data": { - "text/plain": [ - "" + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def build_example(include_redundant_constraint=True):\n", + " m = pyo.ConcreteModel()\n", + "\n", + " m.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", + "\n", + " m.x = pyo.Var(m.I, bounds=(0, 5), initialize=1.0)\n", + "\n", + " m.con1 = pyo.Constraint(expr=m.x[1] + m.x[2] >= 1)\n", + " m.con2 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", + " m.con3 = pyo.Constraint(expr=m.x[2] - 2 * m.x[3] <= 1)\n", + " m.con4 = pyo.Constraint(expr=m.x[1] + m.x[3] >= 1)\n", + " \n", + " if include_redundant_constraint:\n", + " m.con5 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", + "\n", + " m.obj = pyo.Objective(expr=sum(m.x[i] for i in m.I))\n", + " \n", + " return m\n", + "\n", + "m = build_example()" ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dh.check_variable_bounds(tol=1.0)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve the optimization problem" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we can solve the optimization problem. We first set the number of iterations to 50 and then resolve with Ipopt." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Ipopt 3.13.2: max_iter=50\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", - " 1 5.5940309e+00 8.19e+00 2.58e+00 -1.0 4.12e+00 - 3.29e-01 9.99e-02h 1\n", - " 2 4.3446689e+01 2.13e+00 2.10e+00 -1.0 3.67e+00 - 7.64e-01 1.00e+00h 1\n", - " 3 5.0854561e+01 6.80e-01 1.02e+00 -1.0 1.02e+00 0.0 9.34e-01 1.00e+00h 1\n", - " 4 4.9204340e+01 4.73e-01 1.26e+00 -1.0 1.78e+00 - 1.00e+00 1.00e+00f 1\n", - " 5 4.8955163e+01 9.44e-01 1.57e+00 -1.0 4.62e+00 -0.5 1.00e+00 3.51e-01f 2\n", - " 6 4.5850988e+01 8.92e-01 1.56e+00 -1.0 2.14e+00 -0.1 1.00e+00 7.67e-01f 1\n", - " 7 4.5704084e+01 1.19e+00 1.75e+00 -1.0 4.73e+00 - 1.00e+00 2.50e-01f 3\n", - " 8 4.5331671e+01 1.31e+00 3.01e+00 -1.0 1.56e+01 - 7.20e-01 5.98e-02f 4\n", - " 9 4.1480725e+01 7.76e-03 1.98e-01 -1.0 4.94e-01 - 1.00e+00 1.00e+00f 1\n", - "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 10 4.1177584e+01 1.30e-03 8.56e-03 -2.5 1.48e-01 - 1.00e+00 1.00e+00f 1\n", - " 11 4.1157407e+01 7.91e-06 1.09e-05 -3.8 9.88e-03 - 1.00e+00 1.00e+00h 1\n", - " 12 4.1157067e+01 2.59e-10 1.57e-09 -5.7 1.30e-04 - 1.00e+00 1.00e+00h 1\n", - " 13 4.1157063e+01 2.24e-14 1.63e-13 -8.6 1.52e-06 - 1.00e+00 1.00e+00h 1\n", - "\n", - "Number of Iterations....: 13\n", - "\n", - " (scaled) (unscaled)\n", - "Objective...............: 4.1157063321491172e+01 4.1157063321491172e+01\n", - "Dual infeasibility......: 1.6278818285040336e-13 1.6278818285040336e-13\n", - "Constraint violation....: 2.2426505097428162e-14 2.2426505097428162e-14\n", - "Complementarity.........: 2.5065697306197698e-09 2.5065697306197698e-09\n", - "Overall NLP error.......: 2.5065697306197698e-09 2.5065697306197698e-09\n", - "\n", - "\n", - "Number of objective function evaluations = 24\n", - "Number of objective gradient evaluations = 14\n", - "Number of equality constraint evaluations = 24\n", - "Number of inequality constraint evaluations = 24\n", - "Number of equality constraint Jacobian evaluations = 14\n", - "Number of inequality constraint Jacobian evaluations = 14\n", - "Number of Lagrangian Hessian evaluations = 13\n", - "Total CPU secs in IPOPT (w/o function evaluations) = 0.001\n", - "Total CPU secs in NLP function evaluations = 0.000\n", - "\n", - "EXIT: Optimal Solution Found.\n" - ] }, { - "data": { - "text/plain": [ - "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 3, 'Number of variables': 5, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.014306068420410156}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Check for Structural Issues\n", + "\n", + "Before we try to solve the model, the first thing we should do is use the Diagnostics Toolbox to check for any structural issues in our model. The cell below creates an instance of the Diagnostics Toolbox and calls the ``report_structural_issues`` method." ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m, tee=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, Ipopt has no trouble solving this optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check if any constraint residuals are large than 1E-14" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's now inspect the new solution to see which (if any) constraints have residuals larger than 10$^{-14}$." - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " \n", - "All constraints with residuals larger than 1e-14 :\n", - "Count\tName\t|residual|\n", - "0 \t con1 \t 9.97185747309004e-08\n", - "1 \t con2 \t 7.942586144338293e-09\n", - "2 \t con3 \t 2.2426505097428162e-14\n" - ] }, { - "data": { - "text/plain": [ - "dict_keys([, , ])" + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "Model Statistics\n", + "\n", + " Activated Blocks: 1 (Deactivated: 0)\n", + " Free Variables in Activated Constraints: 3 (External: 0)\n", + " Free Variables with only lower bounds: 0\n", + " Free Variables with only upper bounds: 0\n", + " Free Variables with upper and lower bounds: 3\n", + " Fixed Variables in Activated Constraints: 0 (External: 0)\n", + " Activated Equality Constraints: 2 (Deactivated: 0)\n", + " Activated Inequality Constraints: 3 (Deactivated: 0)\n", + " Activated Objectives: 1 (Deactivated: 0)\n", + "\n", + "------------------------------------------------------------------------------------\n", + "2 WARNINGS\n", + "\n", + " WARNING: 1 Degree of Freedom\n", + " WARNING: Structural singularity found\n", + " Under-Constrained Set: 3 variables, 2 constraints\n", + " Over-Constrained Set: 0 variables, 0 constraints\n", + "\n", + "------------------------------------------------------------------------------------\n", + "0 Cautions\n", + "\n", + " No cautions found!\n", + "\n", + "------------------------------------------------------------------------------------\n", + "Suggested next steps:\n", + "\n", + " display_underconstrained_set()\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "dt = DiagnosticsToolbox(m)\n", + "dt.report_structural_issues()" ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dh.check_residuals(tol=1e-14)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, all of the constraints are satisfied, even with this fairly tight tolerance." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify all variables within 1E-5 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, let's check if any of the variables are near their bounds at the new solution." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " \n", - "No variables within 1e-05 (absolute) of their bounds.\n" - ] }, { - "data": { - "text/plain": [ - "" + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, the toolbox is reporting one degree of freedom and a structural singularity. In this case, these are both expected as we are trying to solve an optimization problem with one degree of freedom so we can move on to the next step. \n", + "\n", + "### Evaluate the initial point\n", + "\n", + "Before we try to solve our degenerate model, it can often be useful to check the state of the model at the initial step seen by the solver. The cell below creates an instance of Ipopt and sets the iteration limit to 0 so that it terminates at the initial condition." ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dh.check_variable_bounds(tol=1e-5)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Great, no variables are near their bounds. If a variable was at its bound, it is important the inspect the model and confirm the bound is physically sensible/what you intended." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check the rank of the constraint Jacobian at the solution" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The main feature of Degeneracy Hunter is to check if an optimization problem is poorly formulated. Let's see what happens when we check the rank on a carefully formulated optimization problem:" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Checking rank of Jacobian of equality constraints...\n", - "Model contains 1 equality constraints and 5 variables.\n", - "Only singular value: 2.23606797749979\n" - ] }, { - "data": { - "text/plain": [ - "0" + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.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...............: 3.0000000000000000e+00 3.0000000000000000e+00\n", + "Dual infeasibility......: 1.0000000000000000e+00 1.0000000000000000e+00\n", + "Constraint violation....: 2.0000000000000000e+00 2.0000000000000000e+00\n", + "Complementarity.........: 4.0000000499999997e+00 4.0000000499999997e+00\n", + "Overall NLP error.......: 4.0000000499999997e+00 4.0000000499999997e+00\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" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'warning', 'Message': 'Ipopt 3.13.2\\\\x3a Maximum Number of Iterations Exceeded.', 'Termination condition': 'maxIterations', 'Id': 400, 'Error rc': 0, 'Time': 0.004778385162353516}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver = pyo.SolverFactory(\"ipopt\")\n", + "\n", + "# Specifying an iteration limit of 0 allows us to inspect the initial point\n", + "solver.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", + "solver.solve(m, tee=True)" ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dh.check_rank_equality_constraints()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example 2: Linear Program with Redundant Equality Constraints" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now let's apply Degeneracy Hunter to a poorly formulated optimization problem:\n", - "\n", - "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", - "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", - "& x_1 + x_2 + x_3 = 1 \\\\\n", - "& x_2 - 2 x_3 \\leq 1 \\\\\n", - "& x_1 + x_3 \\geq 1 \\\\\n", - "& x_1 + x_2 + x_3 = 1 \\\\\n", - "\\end{align*} $$\n", - "\n", - "\n", - "Notice the two equality constraints are redundant. This means the constraint qualifications (e.g., LICQ) do not hold which has three important implications:\n", - "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", - "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", - "3. Theoretical convergence properties of optimization algorithms may not hold\n", - "\n", - "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the primary purpose of Degeneracy Hunter. Let's see it in action." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Define the model in Pyomo" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "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 : 3 : {1, 2, 3}\n", - "\n", - "1 Var Declarations\n", - " x : Size=3, Index=I\n", - " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", - " 1 : 0 : 1.0 : 5 : False : False : Reals\n", - " 2 : 0 : 1.0 : 5 : False : False : Reals\n", - " 3 : 0 : 1.0 : 5 : 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[1] + x[2] + x[3]\n", - "\n", - "5 Constraint Declarations\n", - " con1 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : 1.0 : x[1] + x[2] : +Inf : True\n", - " con2 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : 1.0 : x[1] + x[2] + x[3] : 1.0 : True\n", - " con3 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : -Inf : x[2] - 2*x[3] : 1.0 : True\n", - " con4 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : 1.0 : x[1] + x[3] : +Inf : True\n", - " con5 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : 1.0 : x[1] + x[2] + x[3] : 1.0 : True\n", - "\n", - "8 Declarations: I x con1 con2 con3 con4 con5 obj\n" - ] - } - ], - "source": [ - "def example2(with_degenerate_constraint=True):\n", - " \"\"\"Create the Pyomo model for Example 2\n", - "\n", - " Arguments:\n", - " with_degenerate_constraint: Boolean, if True, include the redundant linear constraint\n", - "\n", - " Returns:\n", - " m2: Pyomo model\n", - " \"\"\"\n", - "\n", - " m2 = pyo.ConcreteModel()\n", - "\n", - " m2.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", - "\n", - " m2.x = pyo.Var(m2.I, bounds=(0, 5), initialize=1.0)\n", - "\n", - " m2.con1 = pyo.Constraint(expr=m2.x[1] + m2.x[2] >= 1)\n", - " m2.con2 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", - " m2.con3 = pyo.Constraint(expr=m2.x[2] - 2 * m2.x[3] <= 1)\n", - " m2.con4 = pyo.Constraint(expr=m2.x[1] + m2.x[3] >= 1)\n", - "\n", - " if with_degenerate_constraint:\n", - " m2.con5 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", - "\n", - " m2.obj = pyo.Objective(expr=sum(m2.x[i] for i in m2.I))\n", - "\n", - " m2.pprint()\n", - "\n", - " return m2\n", - "\n", - "\n", - "# Create the Pyomo model for Example 2 including the redundant constraint\n", - "m2 = example2()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Evaluate the initial point" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "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...: 6\n", - "Number of nonzeros in inequality constraint Jacobian.: 6\n", - "Number of nonzeros in Lagrangian Hessian.............: 0\n", - "\n", - "Total number of variables............................: 3\n", - " variables with only lower bounds: 0\n", - " variables with lower and upper bounds: 3\n", - " variables with only upper bounds: 0\n", - "Total number of equality constraints.................: 2\n", - "Total number of inequality constraints...............: 3\n", - " inequality constraints with only lower bounds: 2\n", - " inequality constraints with lower and upper bounds: 0\n", - " inequality constraints with only upper bounds: 1\n", - "\n", - "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 0 3.0000000e+00 2.00e+00 1.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...............: 3.0000000000000000e+00 3.0000000000000000e+00\n", - "Dual infeasibility......: 1.0000000000000000e+00 1.0000000000000000e+00\n", - "Constraint violation....: 2.0000000000000000e+00 2.0000000000000000e+00\n", - "Complementarity.........: 4.0000000499999997e+00 4.0000000499999997e+00\n", - "Overall NLP error.......: 4.0000000499999997e+00 4.0000000499999997e+00\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" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Identify constraints with large residuals at the initial point\n", + "\n", + "With the solution at the initial point, we can then use the Diagnostics Toolbox to check the residuals of all constraints at the initial point." + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "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" - ] + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following constraint(s) have large residuals (>1.0E-05):\n", + "\n", + " con2: 2.00000E+00\n", + " con5: 2.00000E+00\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "# Check for large residuals\n", + "dt.display_constraints_with_large_residuals()" + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "2023-11-02 10:25:10 [WARNING] idaes.core.util.model_diagnostics: DEPRECATED: DegeneracyHunter is being deprecated in favor of the new\n", - "DiagnosticsToolbox. (deprecated in 2.2.0, will be removed in (or\n", - "after) 3.0.0)\n", - "(called from C:\\Users\\dkgun\\AppData\\Local\\Temp\\ipykernel_30328\\1651168461.py:11)\n" - ] - } - ], - "source": [ - "# 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(m2, tee=True)\n", - "\n", - "# Create Degeneracy Hunter object\n", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh2 = DegeneracyHunter(m2, solver=pyo.SolverFactory(\"cbc\"))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify constraints with residuals greater than 0.1 at the initial point" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also check for variables which are close to (or outside of) their bounds." + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - " \n", - "All constraints with residuals larger than 0.1 :\n", - "Count\tName\t|residual|\n", - "0 \t con2 \t 2.0\n", - "1 \t con5 \t 2.0\n" - ] + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following variable(s) have values close to their bounds (abs=1.0E-04, rel=1.0E-04):\n", + "\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "dt.display_variables_near_bounds()" + ] }, { - "data": { - "text/plain": [ - "dict_keys([, ])" + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Solve the optimization problem and extract the solution\n", + "\n", + "There appear to be no significant issues at the initial point, so let us move on to solving the optimization problem." ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "dh2.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve the optimization problem and extract the solution" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now let's solve the optimization problem." - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [ + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Ipopt 3.13.2: max_iter=50\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...: 6\n", - "Number of nonzeros in inequality constraint Jacobian.: 6\n", - "Number of nonzeros in Lagrangian Hessian.............: 0\n", - "\n", - "Total number of variables............................: 3\n", - " variables with only lower bounds: 0\n", - " variables with lower and upper bounds: 3\n", - " variables with only upper bounds: 0\n", - "Total number of equality constraints.................: 2\n", - "Total number of inequality constraints...............: 3\n", - " inequality constraints with only lower bounds: 2\n", - " inequality constraints with lower and upper bounds: 0\n", - " inequality constraints with only upper bounds: 1\n", - "\n", - "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 0 3.0000000e+00 2.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", - " 1 1.3934810e+00 3.93e-01 2.17e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01h 1\n", - " 2 1.0184419e+00 1.84e-02 1.26e-02 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", - " 3 1.0011246e+00 1.12e-03 3.82e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", - " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", - " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", - " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 4.99e-04 - 1.00e+00 5.82e-01h 1\n", - " 7 1.0000788e+00 7.88e-05 4.33e+01 -2.5 1.88e-04 - 1.00e+00 2.94e-01f 2\n", - " 8 1.0000153e+00 1.53e-05 2.19e+01 -2.5 6.85e-05 - 1.00e+00 8.08e-01h 1\n", - " 9 1.0000118e+00 1.18e-05 2.78e+02 -2.5 3.47e-05 - 1.00e+00 2.45e-01f 2\n" - ] + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 2.17e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01h 1\n", + " 2 1.0184419e+00 1.84e-02 1.26e-02 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.82e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 4.99e-04 - 1.00e+00 5.82e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.33e+01 -2.5 1.88e-04 - 1.00e+00 2.94e-01f 2\n", + " 8 1.0000153e+00 1.53e-05 2.19e+01 -2.5 6.85e-05 - 1.00e+00 8.08e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.78e+02 -2.5 3.47e-05 - 1.00e+00 2.45e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000014e+00 1.44e-06 2.83e-08 -2.5 7.79e-06 - 1.00e+00 1.00e+00h 1\n", + " 11 1.0000000e+00 1.39e-09 1.84e-11 -5.7 2.56e-06 - 1.00e+00 1.00e+00h 1\n", + " 12 1.0000000e+00 1.39e-09 1.24e+02 -8.6 1.16e-08 - 1.00e+00 9.54e-07h 21\n", + " 13 9.9999999e-01 9.82e-09 1.14e-13 -8.6 1.33e-08 - 1.00e+00 1.00e+00s 22\n", + " 14 9.9999999e-01 9.96e-09 6.46e-14 -8.6 2.29e-10 - 1.00e+00 1.00e+00s 22\n", + " 15 9.9999999e-01 9.96e-09 1.82e+02 -9.0 4.00e-11 - 1.00e+00 9.54e-07h 21\n", + " 16 9.9999999e-01 1.00e-08 1.52e-13 -9.0 5.89e-11 - 1.00e+00 1.00e+00s 22\n", + "\n", + "Number of Iterations....: 16\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999000238915e-01 9.9999999000238915e-01\n", + "Dual infeasibility......: 1.5188693260016487e-13 1.5188693260016487e-13\n", + "Constraint violation....: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "Complementarity.........: 9.2217032157601989e-10 9.2217032157601989e-10\n", + "Overall NLP error.......: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 111\n", + "Number of objective gradient evaluations = 17\n", + "Number of equality constraint evaluations = 111\n", + "Number of inequality constraint evaluations = 111\n", + "Number of equality constraint Jacobian evaluations = 17\n", + "Number of inequality constraint Jacobian evaluations = 17\n", + "Number of Lagrangian Hessian evaluations = 16\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: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.04887509346008301}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver.options[\"max_iter\"] = 50\n", + "solver.solve(m, tee=True)" + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 10 1.0000014e+00 1.44e-06 2.83e-08 -2.5 7.79e-06 - 1.00e+00 1.00e+00h 1\n", - " 11 1.0000000e+00 1.39e-09 1.84e-11 -5.7 2.56e-06 - 1.00e+00 1.00e+00h 1\n", - " 12 1.0000000e+00 1.39e-09 1.24e+02 -8.6 1.16e-08 - 1.00e+00 9.54e-07h 21\n", - " 13 9.9999999e-01 9.82e-09 1.14e-13 -8.6 1.33e-08 - 1.00e+00 1.00e+00s 22\n", - " 14 9.9999999e-01 9.96e-09 6.46e-14 -8.6 2.29e-10 - 1.00e+00 1.00e+00s 22\n", - " 15 9.9999999e-01 9.96e-09 1.82e+02 -9.0 4.00e-11 - 1.00e+00 9.54e-07h 21\n", - " 16 9.9999999e-01 1.00e-08 1.52e-13 -9.0 5.89e-11 - 1.00e+00 1.00e+00s 22\n", - "\n", - "Number of Iterations....: 16\n", - "\n", - " (scaled) (unscaled)\n", - "Objective...............: 9.9999999000238915e-01 9.9999999000238915e-01\n", - "Dual infeasibility......: 1.5188693260016487e-13 1.5188693260016487e-13\n", - "Constraint violation....: 9.9976108502985994e-09 9.9976108502985994e-09\n", - "Complementarity.........: 9.2217032157601989e-10 9.2217032157601989e-10\n", - "Overall NLP error.......: 9.9976108502985994e-09 9.9976108502985994e-09\n", - "\n", - "\n", - "Number of objective function evaluations = 111\n", - "Number of objective gradient evaluations = 17\n", - "Number of equality constraint evaluations = 111\n", - "Number of inequality constraint evaluations = 111\n", - "Number of equality constraint Jacobian evaluations = 17\n", - "Number of inequality constraint Jacobian evaluations = 17\n", - "Number of Lagrangian Hessian evaluations = 16\n", - "Total CPU secs in IPOPT (w/o function evaluations) = 0.003\n", - "Total CPU secs in NLP function evaluations = 0.001\n", - "\n", - "EXIT: Optimal Solution Found.\n" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of line search evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution.\n", + "\n", + "The cell below shows the values of the variables at the solution." + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "x[ 1 ]= 1.0000000099975996\n", - "x[ 2 ]= 0.0\n", - "x[ 3 ]= 0.0\n" - ] - } - ], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2, tee=True)\n", - "\n", - "for i in m2.I:\n", - " print(\"x[\", i, \"]=\", m2.x[i]())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of linesearch evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check the rank of the Jacobian of the equality constraints" - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [ + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000099975996 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 0.0 : 5 : False : False : Reals\n" + ] + } + ], + "source": [ + "m.x.display()" + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Checking rank of Jacobian of equality constraints...\n", - "Model contains 2 equality constraints and 3 variables.\n", - "Computing the 1 smallest singular value(s)\n", - "Smallest singular value(s):\n", - "0.000E+00\n" - ] - } - ], - "source": [ - "n_deficient = dh2.check_rank_equality_constraints()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "A singular value near 0 indicates the Jacobian of the equality constraints is rank deficient. For each near-zero singular value, there is likely one degenerate constraint." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify candidate degenerate constraints" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Degeneracy Hunter first identifies candidate degenerate equations." - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Check the rank of the Jacobian of the equality constraints\n", + "\n", + "One way to check for the presence of degenerate equations is to calculate the rank of the Jacobian matrix of the equality constraints. A singular value near 0 indicates the Jacobian is rank deficient, and for each near-zero singular value there is likely one degenerate constraint.\n", + "\n", + "The Diagnostics Toolbox contains a Singular Value Decomposition (SVD) toolbox which can be used to determine the number of small singular values in the Jacobian for a model as shown below." + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "*** Searching for a Single Degenerate Set ***\n", - "Building MILP model...\n", - "2 Set Declarations\n", - " C : Size=1, Index=None, Ordered=Insertion\n", - " Key : Dimen : Domain : Size : Members\n", - " None : 1 : Any : 2 : {0, 1}\n", - " V : Size=1, Index=None, Ordered=Insertion\n", - " Key : Dimen : Domain : Size : Members\n", - " None : 1 : Any : 3 : {0, 1, 2}\n", - "\n", - "4 Var Declarations\n", - " abs_nu : Size=2, Index=C\n", - " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", - " 0 : 0 : None : 100000.00001 : False : True : Reals\n", - " 1 : 0 : None : 100000.00001 : False : True : Reals\n", - " nu : Size=2, Index=C\n", - " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", - " 0 : -100000.00001 : 1.0 : 100000.00001 : False : False : Reals\n", - " 1 : -100000.00001 : 1.0 : 100000.00001 : False : False : Reals\n", - " y_neg : Size=2, Index=C\n", - " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", - " 0 : 0 : None : 1 : False : True : Binary\n", - " 1 : 0 : None : 1 : False : True : Binary\n", - " y_pos : Size=2, Index=C\n", - " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", - " 0 : 0 : None : 1 : False : True : Binary\n", - " 1 : 0 : None : 1 : False : True : Binary\n", - "\n", - "1 Constraint Declarations\n", - " pos_xor_neg : Size=0, Index=C, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - "\n", - "7 Declarations: C V nu y_pos y_neg abs_nu pos_xor_neg\n", - "Solving MILP model...\n" - ] + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "\n", + "Number of Singular Values less than 1.0E-06 is 1\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "svd = dt.prepare_svd_toolbox()\n", + "svd.display_rank_of_equality_constraints()" + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Welcome to the CBC MILP Solver \n", - "Version: 2.10.10 \n", - "Build Date: Jun 8 2023 \n", - "\n", - "command line - C:\\Users\\dkgun\\AppData\\Local\\idaes\\bin\\cbc.exe -printingOptions all -import C:\\Users\\dkgun\\AppData\\Local\\Temp\\tmp1neqsh_k.pyomo.lp -stat=1 -solve -solu C:\\Users\\dkgun\\AppData\\Local\\Temp\\tmp1neqsh_k.pyomo.soln (default strategy 1)\n", - "Option for printingOptions changed from normal to all\n", - "Presolve 9 (-7) rows, 7 (-1) columns and 20 (-14) elements\n", - "Statistics for presolved model\n", - "Original problem has 4 integers (4 of which binary)\n", - "Presolved problem has 4 integers (4 of which binary)\n", - "==== 5 zero objective 2 different\n", - "5 variables have objective of 0\n", - "2 variables have objective of 1\n", - "==== absolute objective values 2 different\n", - "5 variables have objective of 0\n", - "2 variables have objective of 1\n", - "==== for integers 4 zero objective 1 different\n", - "4 variables have objective of 0\n", - "==== for integers absolute objective values 1 different\n", - "4 variables have objective of 0\n", - "===== end objective counts\n", - "\n", - "\n", - "Problem has 9 rows, 7 columns (2 with objective) and 20 elements\n", - "Column breakdown:\n", - "0 of type 0.0->inf, 2 of type 0.0->up, 0 of type lo->inf, \n", - "1 of type lo->up, 0 of type free, 0 of type fixed, \n", - "0 of type -inf->0.0, 0 of type -inf->up, 4 of type 0.0->1.0 \n", - "Row breakdown:\n", - "0 of type E 0.0, 0 of type E 1.0, 0 of type E -1.0, \n", - "0 of type E other, 0 of type G 0.0, 1 of type G 1.0, \n", - "0 of type G other, 4 of type L 0.0, 0 of type L 1.0, \n", - "4 of type L other, 0 of type Range 0.0->1.0, 0 of type Range other, \n", - "0 of type Free \n", - "Continuous objective value is 0 - 0.00 seconds\n", - "Cgl0003I 0 fixed, 0 tightened bounds, 2 strengthened rows, 0 substitutions\n", - "Cgl0004I processed model has 9 rows, 7 columns (4 integer (4 of which binary)) and 22 elements\n", - "Cbc0038I Initial state - 2 integers unsatisfied sum - 1\n", - "Cbc0038I Pass 1: suminf. 0.00000 (0) obj. 1.99999e-05 iterations 3\n", - "Cbc0038I Solution found of 1.99999e-05\n", - "Cbc0038I Relaxing continuous gives 2e-05\n", - "Cbc0038I Rounding solution of 1.99999e-05 is better than previous of 2e-05\n", - "\n", - "Cbc0038I Before mini branch and bound, 2 integers at bound fixed and 0 continuous\n", - "Cbc0038I Full problem 9 rows 7 columns, reduced to 7 rows 5 columns\n", - "Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)\n", - "Cbc0038I Round again with cutoff of 8.99991e-06\n", - "Cbc0038I Pass 2: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 2\n", - "Cbc0038I Pass 3: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 4: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 1\n", - "Cbc0038I Pass 5: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 6: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 7: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 8: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 9: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 10: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 11: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 12: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 1\n", - "Cbc0038I Pass 13: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 14: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 1\n", - "Cbc0038I Pass 15: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 16: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 17: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 18: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 5\n", - "Cbc0038I Pass 19: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 1\n", - "Cbc0038I Pass 20: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 21: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 22: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 23: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 24: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 25: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 26: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 2\n", - "Cbc0038I Pass 27: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 2\n", - "Cbc0038I Pass 28: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 0\n", - "Cbc0038I Pass 29: suminf. 1.00000 (3) obj. 8.99991e-06 iterations 3\n", - "Cbc0038I Pass 30: suminf. 0.55000 (2) obj. 8.99991e-06 iterations 2\n", - "Cbc0038I Pass 31: suminf. 1.00000 (3) obj. 8.99991e-06 iterations 2\n", - "Cbc0038I Rounding solution of 8.99991e-06 is better than previous of 1.99999e-05\n", - "\n", - "Cbc0038I Before mini branch and bound, 0 integers at bound fixed and 0 continuous\n", - "Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)\n", - "Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 8.99991e-06 - took 0.00 seconds\n", - "Cbc0012I Integer solution of 2e-05 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)\n", - "Cbc0038I Full problem 9 rows 7 columns, reduced to 8 rows 6 columns\n", - "Cbc0006I The LP relaxation is infeasible or too expensive\n", - "Cbc0031I 2 added rows had average density of 2\n", - "Cbc0013I At root node, 2 cuts changed objective from 0 to 1e-05 in 3 passes\n", - "Cbc0014I Cut generator 0 (Probing) - 2 row cuts average 1.0 elements, 2 column cuts (2 active) in 0.000 seconds - new frequency is 1\n", - "Cbc0014I Cut generator 1 (Gomory) - 2 row cuts average 2.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1\n", - "Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100\n", - "Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100\n", - "Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100\n", - "Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100\n", - "Cbc0014I Cut generator 6 (TwoMirCuts) - 2 row cuts average 2.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1\n", - "Cbc0001I Search completed - best objective 2.000000000000001e-05, took 6 iterations and 0 nodes (0.00 seconds)\n", - "Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost\n", - "Cuts at root node changed objective from 0 to 1e-05\n", - "Probing was tried 3 times and created 4 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "Gomory was tried 2 times and created 2 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "Knapsack was tried 2 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "Clique was tried 2 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "MixedIntegerRounding2 was tried 2 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "FlowCover was tried 2 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "TwoMirCuts was tried 2 times and created 2 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "ZeroHalf was tried 1 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)\n", - "\n", - "Result - Optimal solution found\n", - "\n", - "Objective value: 0.00002000\n", - "Enumerated nodes: 0\n", - "Total iterations: 6\n", - "Time (CPU seconds): 0.01\n", - "Time (Wallclock seconds): 0.00\n", - "\n", - "Total time (CPU seconds): 0.01 (Wallclock seconds): 0.01\n", - "\n" - ] - } - ], - "source": [ - "ds2 = dh2.find_candidate_equations(verbose=True, tee=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Find irreducible degenerate sets (IDS)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Next, Degeneracy Hunter enumerates through the candidate equations. For each candidate equation, Degenerate Hunter solves a MILP to compute the irreducible degenerate set that must contain the candidate equation." - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As can be seen above, there is one singular value less than 1e-6, suggesting we have one degenerate constraint." + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "*** Searching for Irreducible Degenerate Sets ***\n", - "Building MILP model...\n", - "Solving MILP 1 of 2 ...\n" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Finding Irreducible Degenerate Sets (IDS)\n", + "\n", + "Next, we can use Degeneracy Hunter to identify irreducible degenerate sets; that is, the smallest set of constraints that contain a redundant constraint.\n", + "\n", + "Degeneracy Hunter first identifies candidate degenerate equations by solving a mixed integer linear program (MILP). It then iterates through the candidate equations and solves a second MILP for each to compute the irreducible degenerate set that must contain the candidate equation.\n", + "\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver, users can specify their solver of choice via a configuration argument. The default option is SCIP (https://scipopt.org/) which is an open-source solver.\n", + "
\n", + "\n", + "The cell below shows how to create an instance of Degeneracy Hunter and generate a report of the irreducible degenerate sets." + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Solving MILP 2 of 2 ...\n" - ] + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Candidate Equations\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving Candidates MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Irreducible Degenerate Sets\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model to compute irreducible degenerate set.\n", + "Solving MILP 1 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con2.\n", + "Solving MILP 2 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con5.\n", + "====================================================================================\n", + "Irreducible Degenerate Sets\n", + "\n", + " Irreducible Degenerate Set 0\n", + " nu Constraint Name\n", + " 1.0 con2\n", + " -1.0 con5\n", + "\n", + " Irreducible Degenerate Set 1\n", + " nu Constraint Name\n", + " -1.0 con2\n", + " 1.0 con5\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "dh = dt.prepare_degeneracy_hunter()\n", + "dh.report_irreducible_degenerate_sets()" + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n", - "Irreducible Degenerate Set" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As can be seen above, Degeneracy Hunter identified two constraints as potentially redundant (the duplicate constraints ``con2`` and ``con5``). Degeneracy Hunter then reports an Irreducible Degenerate Set for each of these constraints, which in the case as identical. More complex models may have partially overlapping IDS's.\n", + "\n", + "These results show us that ``con2`` and ``con5`` are mutually redundant; i.e., we can eliminate one of these with no effect on the model solution." + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - " 0\n", - "nu\tConstraint Name\n", - "1.0 \t con2\n", - "-1.0 \t con5\n", - "\n", - "Irreducible Degenerate Set 1\n", - "nu\tConstraint Name\n", - "-1.0 \t con2\n", - "1.0 \t con5\n" - ] - } - ], - "source": [ - "ids = dh2.find_irreducible_degenerate_sets(verbose=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Reformulate Example 2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now let's reformulate the model by skipping/removing the redundant equality constraint:\n", - "\n", - "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", - "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", - "& x_1 + x_2 + x_3 = 1 \\\\\n", - "& x_2 - 2 x_3 \\leq 1 \\\\\n", - "& x_1 + x_3 \\geq 1\n", - "\\end{align*} $$" - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "metadata": {}, - "outputs": [ + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Reformulate Model to Remove Redundant Constraint\n", + "\n", + "Now let's reformulate the model by removing the redundant equality constraint ``con5`` and solve the reformulated model." + ] + }, { - "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 : 3 : {1, 2, 3}\n", - "\n", - "1 Var Declarations\n", - " x : Size=3, Index=I\n", - " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", - " 1 : 0 : 1.0 : 5 : False : False : Reals\n", - " 2 : 0 : 1.0 : 5 : False : False : Reals\n", - " 3 : 0 : 1.0 : 5 : 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[1] + x[2] + x[3]\n", - "\n", - "4 Constraint Declarations\n", - " con1 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : 1.0 : x[1] + x[2] : +Inf : True\n", - " con2 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : 1.0 : x[1] + x[2] + x[3] : 1.0 : True\n", - " con3 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : -Inf : x[2] - 2*x[3] : 1.0 : True\n", - " con4 : Size=1, Index=None, Active=True\n", - " Key : Lower : Body : Upper : Active\n", - " None : 1.0 : x[1] + x[3] : +Inf : True\n", - "\n", - "7 Declarations: I x con1 con2 con3 con4 obj\n" - ] - } - ], - "source": [ - "m2b = example2(with_degenerate_constraint=False)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve the reformulated model" - ] - }, - { - "cell_type": "code", - "execution_count": 19, - "metadata": {}, - "outputs": [ + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 1\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 6.30e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 1.34e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01f 1\n", + " 2 1.0184419e+00 1.84e-02 9.15e-03 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.88e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 5.00e-04 - 1.00e+00 5.81e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.34e+01 -2.5 1.88e-04 - 1.00e+00 2.93e-01f 2\n", + " 8 1.0000154e+00 1.54e-05 2.26e+01 -2.5 6.89e-05 - 1.00e+00 8.04e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.98e+02 -2.5 3.64e-05 - 1.00e+00 2.33e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000020e+00 2.04e-06 1.33e+02 -2.5 9.72e-06 - 1.00e+00 8.28e-01h 1\n", + " 11 1.0000016e+00 1.61e-06 2.26e+03 -2.5 4.79e-06 - 1.00e+00 2.12e-01f 2\n", + " 12 1.0000002e+00 2.46e-07 8.72e+02 -2.5 1.20e-06 - 1.00e+00 8.47e-01h 1\n", + " 13 1.0000002e+00 1.95e-07 1.71e+04 -2.5 7.02e-07 - 1.00e+00 2.09e-01f 2\n", + " 14 1.0000000e+00 1.54e-08 3.23e+03 -2.5 1.50e-07 - 1.00e+00 9.21e-01h 1\n", + " 15 1.0000000e+00 1.15e-08 9.99e+04 -2.5 6.89e-08 - 1.00e+00 2.54e-01f 2\n", + " 16 1.0000000e+00 2.22e-16 2.83e-08 -2.5 8.21e-09 - 1.00e+00 1.00e+00h 1\n", + " 17 1.0000000e+00 2.22e-16 4.14e-11 -8.6 8.25e-16 - 1.00e+00 1.00e+00 0\n", + "\n", + "Number of Iterations....: 17\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999999999978e-01 9.9999999999999978e-01\n", + "Dual infeasibility......: 4.1425156707686206e-11 4.1425156707686206e-11\n", + "Constraint violation....: 2.2204460492503131e-16 2.2204460492503131e-16\n", + "Complementarity.........: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "Overall NLP error.......: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 23\n", + "Number of objective gradient evaluations = 18\n", + "Number of equality constraint evaluations = 23\n", + "Number of inequality constraint evaluations = 23\n", + "Number of equality constraint Jacobian evaluations = 18\n", + "Number of inequality constraint Jacobian evaluations = 18\n", + "Number of Lagrangian Hessian evaluations = 17\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.002\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 4, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.006506204605102539}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m2 = build_example(include_redundant_constraint=False)\n", + "\n", + "solver.solve(m2, tee=True)" + ] + }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "Ipopt 3.13.2: max_iter=50\n" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let us check to see if the optimal solution has changed." + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "\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.............: 0\n", - "\n", - "Total number of variables............................: 3\n", - " variables with only lower bounds: 0\n", - " variables with lower and upper bounds: 3\n", - " variables with only upper bounds: 0\n", - "Total number of equality constraints.................: 1\n", - "Total number of inequality constraints...............: 3\n", - " inequality constraints with only lower bounds: 2\n", - " inequality constraints with lower and upper bounds: 0\n", - " inequality constraints with only upper bounds: 1\n", - "\n", - "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 0 3.0000000e+00 2.00e+00 6.30e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", - " 1 1.3934810e+00 3.93e-01 1.34e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01f 1\n", - " 2 1.0184419e+00 1.84e-02 9.15e-03 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", - " 3 1.0011246e+00 1.12e-03 3.88e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", - " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", - " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", - " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 5.00e-04 - 1.00e+00 5.81e-01h 1\n", - " 7 1.0000788e+00 7.88e-05 4.34e+01 -2.5 1.88e-04 - 1.00e+00 2.93e-01f 2\n", - " 8 1.0000154e+00 1.54e-05 2.26e+01 -2.5 6.89e-05 - 1.00e+00 8.04e-01h 1\n", - " 9 1.0000118e+00 1.18e-05 2.98e+02 -2.5 3.64e-05 - 1.00e+00 2.33e-01f 2\n", - "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", - " 10 1.0000020e+00 2.04e-06 1.33e+02 -2.5 9.72e-06 - 1.00e+00 8.28e-01h 1\n", - " 11 1.0000016e+00 1.61e-06 2.26e+03 -2.5 4.79e-06 - 1.00e+00 2.12e-01f 2\n", - " 12 1.0000002e+00 2.46e-07 8.72e+02 -2.5 1.20e-06 - 1.00e+00 8.47e-01h 1\n", - " 13 1.0000002e+00 1.95e-07 1.71e+04 -2.5 7.02e-07 - 1.00e+00 2.09e-01f 2\n", - " 14 1.0000000e+00 1.54e-08 3.23e+03 -2.5 1.50e-07 - 1.00e+00 9.21e-01h 1\n", - " 15 1.0000000e+00 1.15e-08 9.99e+04 -2.5 6.89e-08 - 1.00e+00 2.54e-01f 2\n", - " 16 1.0000000e+00 2.22e-16 2.83e-08 -2.5 8.21e-09 - 1.00e+00 1.00e+00h 1\n", - " 17 1.0000000e+00 2.22e-16 4.14e-11 -8.6 8.25e-16 - 1.00e+00 1.00e+00 0\n", - "\n", - "Number of Iterations....: 17\n", - "\n", - " (scaled) (unscaled)\n", - "Objective...............: 9.9999999999999978e-01 9.9999999999999978e-01\n", - "Dual infeasibility......: 4.1425156707686206e-11 4.1425156707686206e-11\n", - "Constraint violation....: 2.2204460492503131e-16 2.2204460492503131e-16\n", - "Complementarity.........: 2.6658012082875325e-09 2.6658012082875325e-09\n", - "Overall NLP error.......: 2.6658012082875325e-09 2.6658012082875325e-09\n", - "\n", - "\n", - "Number of objective function evaluations = 23\n", - "Number of objective gradient evaluations = 18\n", - "Number of equality constraint evaluations = 23\n", - "Number of inequality constraint evaluations = 23\n", - "Number of equality constraint Jacobian evaluations = 18\n", - "Number of inequality constraint Jacobian evaluations = 18\n", - "Number of Lagrangian Hessian evaluations = 17\n", - "Total CPU secs in IPOPT (w/o function evaluations) = 0.002\n", - "Total CPU secs in NLP function evaluations = 0.000\n", - "\n", - "EXIT: Optimal Solution Found.\n" - ] + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000005924994 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 4.55844318120227e-11 : 5 : False : False : Reals\n" + ] + } + ], + "source": [ + "m2.x.display()" + ] }, { - "name": "stdout", - "output_type": "stream", - "text": [ - "x[ 1 ]= 5.061595738300216\n", - "x[ 2 ]= 9.184026563775058e-26\n", - "x[ 3 ]= -3.4867300513803765\n" - ] + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We get the same answer as before, but careful inspection of the Ipopt output reveals a subtle improvement. Notice `ls` is only 1 or 2 for all of the iterations, in contrast to more than 20 for the original model. This means Ipopt is taking (nearly) full steps for all iterations.\n", + "\n", + "Let's also compare the number of function evaluations.\n", + "\n", + "Original model (using Ipopt 3.13.2 with `ma27`):\n", + "```\n", + "Number of objective function evaluations = 111\n", + "Number of objective gradient evaluations = 17\n", + "Number of equality constraint evaluations = 111\n", + "Number of inequality constraint evaluations = 111\n", + "Number of equality constraint Jacobian evaluations = 17\n", + "Number of inequality constraint Jacobian evaluations = 17\n", + "Number of Lagrangian Hessian evaluations = 16\n", + "```\n", + "\n", + "Reformulated model (using Ipopt 3.13.2 with `ma27`):\n", + "```\n", + "Number of objective function evaluations = 23\n", + "Number of objective gradient evaluations = 18\n", + "Number of equality constraint evaluations = 23\n", + "Number of inequality constraint evaluations = 23\n", + "Number of equality constraint Jacobian evaluations = 18\n", + "Number of inequality constraint Jacobian evaluations = 18\n", + "Number of Lagrangian Hessian evaluations = 17\n", + "```\n", + "\n", + "Removing a **single redundant constraint** reduced the number of objective and constraint evaluations by a **factor of 5**!\n", + "\n", + "Often degenerate equations have a much worse impact on large-scale problems; for example, degenerate equations can cause Ipopt to require many more iterations or terminate at an infeasible point.\n" + ] + } + ], + "metadata": { + "celltoolbar": "Tags", + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.5" } - ], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2b, tee=True)\n", - "\n", - "for i in m2b.I:\n", - " print(\"x[\", i, \"]=\", m.x[i]())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We get the same answer as before, but careful inspection of the Ipopt output reveals a subtle improvement. Notice `ls` is only 1 or 2 for all of the iterations, in contrast to more than 20 for the original model. This means Ipopt is taking (nearly) full steps for all iterations.\n", - "\n", - "Let's also compare the number of function evaluations.\n", - "\n", - "Original model (using Ipopt 3.13.2 with `ma27`):\n", - "```\n", - "Number of objective function evaluations = 111\n", - "Number of objective gradient evaluations = 17\n", - "Number of equality constraint evaluations = 111\n", - "Number of inequality constraint evaluations = 111\n", - "Number of equality constraint Jacobian evaluations = 17\n", - "Number of inequality constraint Jacobian evaluations = 17\n", - "Number of Lagrangian Hessian evaluations = 16\n", - "```\n", - "\n", - "Reformulated model (using Ipopt 3.13.2 with `ma27`):\n", - "```\n", - "Number of objective function evaluations = 23\n", - "Number of objective gradient evaluations = 18\n", - "Number of equality constraint evaluations = 23\n", - "Number of inequality constraint evaluations = 23\n", - "Number of equality constraint Jacobian evaluations = 18\n", - "Number of inequality constraint Jacobian evaluations = 18\n", - "Number of Lagrangian Hessian evaluations = 17\n", - "```\n", - "\n", - "Removing a **single redundant constraint** reduced the number of objective and constraint evaluations by a **factor of 5**!\n", - "\n", - "Often degenerate equations have a much worse impact on large-scale problems; for example, degenerate equations can cause Ipopt to require many more iterations or terminate at an infeasible point.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "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.11.5" - } - }, - "nbformat": 4, - "nbformat_minor": 3 + "nbformat": 4, + "nbformat_minor": 3 } \ No newline at end of file diff --git a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_test.ipynb b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_test.ipynb index 1b4e368e..569b52c5 100644 --- a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_test.ipynb +++ b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": { "tags": [ "header", @@ -34,324 +34,97 @@ "Maintainer: Andrew Lee \n", "Updated: 2023-06-01 \n", "\n", - "This notebook shows how to use the following Degeneracy Hunter features using two motivating examples:\n", + "> Alexander W. Dowling and Lorenz T. Biegler (2015). Degeneracy Hunter: An Algorithm for Determining Irreducible Sets of Degenerate Constraints in Mathematical Programs. 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering. Ed. by Krist V. Gernaey, Jakob K. Huusom, and Raqul Gani. Computer-Aided Chemical Engineering, Vol. 37, p. 809 \u2013 814. [link to PDF](https://dowlinglab.nd.edu/assets/447116/2015_degeneracy_hunter_an_algorithm_for_determining_irreducible_sets_of_degenerate_constraints_in_mathematical_prog.pdf)\n", + "\n", + "This notebook shows how to use the Degeneracy Hunter tool which is part of the IDAES Diagnostics Toolbox features using two motivating examples:\n", + "\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": "markdown", - "metadata": {}, - "source": [ - "## Setup\n", + "* Demonstrates the performance benefits in Ipopt from removing a single redundant constraint\n", "\n", - "We start by importing Pyomo and Degeneracy Hunter." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import pyomo.environ as pyo\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver. This example uses the open-source solver SCIP (https://scipopt.org/), however users may specify another solver of their choice if they have one available.\n", + "
\n", "\n", - "from idaes.core.util.model_diagnostics import DegeneracyHunter" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example 1: Well-Behaved Nonlinear Program" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Consider the following \"well-behaved\" nonlinear optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "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": {}, - "source": [ - "### Define the model in Pyomo" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We start by defining the optimization problem in Pyomo." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "m = pyo.ConcreteModel()\n", + "## What does Degeneracy Mean?\n", "\n", - "m.I = pyo.Set(initialize=[i for i in range(5)])\n", + "In simple terms, a degenerate model contains one or more constraints that do not add any additional information that is not already included in other constraints in the model. The simplest example of this is the case where a constraint in the model is duplicated. For example, consider the model below:\n", "\n", - "m.x = pyo.Var(m.I, bounds=(-10, 10), initialize=1.0)\n", + "$$\\begin{align*}\n", + "& x_1 + x_2 = 1 \\\\\n", + "& x_1 + x_2 = 1 \\\\\n", + "\\end{align*} $$\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", + "Notice the two equality constraints are identical and thus redundant. This means the constraint qualifications that solvers rely on (e.g., [Linear Independence Constraint Qualification or LICQ](https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions#Regularity_conditions_(or_constraint_qualifications))) do not hold which has three important implications:\n", "\n", - "m.obj = pyo.Objective(expr=sum(m.x[i] ** 2 for i in m.I))\n", + "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", + "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", + "3. Theoretical convergence properties of optimization algorithms may not hold\n", "\n", - "m.pprint()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Evaluate the initial point" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": 3, - "metadata": {}, - "outputs": [], - "source": [ - "# Specify Ipopt as the solver\n", - "opt = pyo.SolverFactory(\"ipopt\")\n", + "Put another way, for a deterministic (square) problem we required that the number of equality constraints be equal to the number of free variables. In this case, there appear to be two equality constraints but in reality there is effectively only one as the two constraints are the same. Thus, the problem is not well defined and there exist a family of solutions that can satisfy the constraints given; any combination of values for ``x1`` and ``x2`` that sum to one will satisfy the constraints as written.\n", "\n", - "# Specifying an iteration limit of 0 allows us to inspect the initial point\n", - "opt.options[\"max_iter\"] = 0\n", + "In practice, degeneracies are rarely as straight-forward as a duplicated constraint. A more common occurrence is a set of linearly dependent constraints; that is a system of constraints where one constraint in the set can be formed by a linear combination of other constraints in the set. For example:\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", + "$$\\begin{align*}\n", + "& x_1 + x_2 + x_3 = 1 \\\\\n", + "& 2x_1 + 2x_2 + 2x_3 = 2 \\\\\n", + "\\end{align*} $$\n", "\n", - "# Create Degeneracy Hunter object\n", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh = DegeneracyHunter(m, solver=pyo.SolverFactory(\"cbc\"))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "source": [ - "### Identify the constraint residuals larger than 0.1" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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", + "Here, the second constraint can be formed by multiplying the first constraint by two. Another common example of linearly dependent constraints is the combination of a set of material balances for all components in a system and an overall material balance. In this case, the overall material balances can be formed by adding all the individual component balances together. Depending on the number of components being modeled, this system of constraints can be quite large, but the end result is the same; the model effectively has one less constraint than it would appear to have.\n", "\n", - "The following line of code will print out all constraints with residuals larger than `0.1`:" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Important: Ipopt does several preprocessing steps when we executed it with zero iterations. When checking the initial point, it is strongly recommended to call Ipopt with zero iterations first. Otherwise, you will not be analyzing the initial point Ipopt starts with." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify all variables within 1 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Another common question when developing optimization models is, \"Which variables are within their bounds by a given tolerance?\" Below is the syntax:" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_variable_bounds(tol=1.0)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve the optimization problem" + "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the purpose of Degeneracy Hunter. Let's see it in action." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now we can solve the optimization problem. We first set the number of iterations to 50 and then resolve with Ipopt." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m, tee=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, Ipopt has no trouble solving this optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check if any constraint residuals are large than 1E-14" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's now inspect the new solution to see which (if any) constraints have residuals larger than 10$^{-14}$." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_residuals(tol=1e-14)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, all of the constraints are satisfied, even with this fairly tight tolerance." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify all variables within 1E-5 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, let's check if any of the variables are near their bounds at the new solution." + "## Setup\n", + "\n", + "We start by importing Pyomo and the IDAES Diagnostics Toolbox." ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "dh.check_variable_bounds(tol=1e-5)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Great, no variables are near their bounds. If a variable was at its bound, it is important the inspect the model and confirm the bound is physically sensible/what you intended." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check the rank of the constraint Jacobian at the solution" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The main feature of Degeneracy Hunter is to check if an optimization problem is poorly formulated. Let's see what happens when we check the rank on a carefully formulated optimization problem:" + "import pyomo.environ as pyo\n", + "\n", + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox" ] }, { "cell_type": "code", - "execution_count": 9, - "metadata": {}, + "execution_count": 3, + "metadata": { + "tags": [ + "testing" + ] + }, "outputs": [], "source": [ - "dh.check_rank_equality_constraints()" + "from idaes.core.util.testing import _enable_scip_solver_for_testing\n", + "scip_solver = pyo.SolverFactory(\"scip\")\n", + "if not scip_solver.available():\n", + " _enable_scip_solver_for_testing()\n", + "assert scip_solver.available()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Example 2: Linear Program with Redundant Equality Constraints" + "## Example: Linear Program with Redundant Equality Constraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's apply Degeneracy Hunter to a poorly formulated optimization problem:\n", + "Let's apply these tools to a poorly formulated optimization problem:\n", "\n", "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", @@ -361,256 +134,661 @@ "& x_1 + x_2 + x_3 = 1 \\\\\n", "\\end{align*} $$\n", "\n", + "As you can see, this is similar to our example above with duplicated equality constraints.\n", "\n", - "Notice the two equality constraints are redundant. This means the constraint qualifications (e.g., LICQ) do not hold which has three important implications:\n", - "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", - "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", - "3. Theoretical convergence properties of optimization algorithms may not hold\n", - "\n", - "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the primary purpose of Degeneracy Hunter. Let's see it in action." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Define the model in Pyomo" + "The cell below shows how to construct this model in Pyomo." ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ - "def example2(with_degenerate_constraint=True):\n", - " \"\"\"Create the Pyomo model for Example 2\n", - "\n", - " Arguments:\n", - " with_degenerate_constraint: Boolean, if True, include the redundant linear constraint\n", - "\n", - " Returns:\n", - " m2: Pyomo model\n", - " \"\"\"\n", - "\n", - " m2 = pyo.ConcreteModel()\n", - "\n", - " m2.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", - "\n", - " m2.x = pyo.Var(m2.I, bounds=(0, 5), initialize=1.0)\n", - "\n", - " m2.con1 = pyo.Constraint(expr=m2.x[1] + m2.x[2] >= 1)\n", - " m2.con2 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", - " m2.con3 = pyo.Constraint(expr=m2.x[2] - 2 * m2.x[3] <= 1)\n", - " m2.con4 = pyo.Constraint(expr=m2.x[1] + m2.x[3] >= 1)\n", - "\n", - " if with_degenerate_constraint:\n", - " m2.con5 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", + "def build_example(include_redundant_constraint=True):\n", + " m = pyo.ConcreteModel()\n", "\n", - " m2.obj = pyo.Objective(expr=sum(m2.x[i] for i in m2.I))\n", + " m.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", "\n", - " m2.pprint()\n", + " m.x = pyo.Var(m.I, bounds=(0, 5), initialize=1.0)\n", "\n", - " return m2\n", + " m.con1 = pyo.Constraint(expr=m.x[1] + m.x[2] >= 1)\n", + " m.con2 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", + " m.con3 = pyo.Constraint(expr=m.x[2] - 2 * m.x[3] <= 1)\n", + " m.con4 = pyo.Constraint(expr=m.x[1] + m.x[3] >= 1)\n", + " \n", + " if include_redundant_constraint:\n", + " m.con5 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", "\n", + " m.obj = pyo.Objective(expr=sum(m.x[i] for i in m.I))\n", + " \n", + " return m\n", "\n", - "# Create the Pyomo model for Example 2 including the redundant constraint\n", - "m2 = example2()" + "m = build_example()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Evaluate the initial point" + "### Check for Structural Issues\n", + "\n", + "Before we try to solve the model, the first thing we should do is use the Diagnostics Toolbox to check for any structural issues in our model. The cell below creates an instance of the Diagnostics Toolbox and calls the ``report_structural_issues`` method." ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 5, "metadata": {}, - "outputs": [], - "source": [ + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "Model Statistics\n", + "\n", + " Activated Blocks: 1 (Deactivated: 0)\n", + " Free Variables in Activated Constraints: 3 (External: 0)\n", + " Free Variables with only lower bounds: 0\n", + " Free Variables with only upper bounds: 0\n", + " Free Variables with upper and lower bounds: 3\n", + " Fixed Variables in Activated Constraints: 0 (External: 0)\n", + " Activated Equality Constraints: 2 (Deactivated: 0)\n", + " Activated Inequality Constraints: 3 (Deactivated: 0)\n", + " Activated Objectives: 1 (Deactivated: 0)\n", + "\n", + "------------------------------------------------------------------------------------\n", + "2 WARNINGS\n", + "\n", + " WARNING: 1 Degree of Freedom\n", + " WARNING: Structural singularity found\n", + " Under-Constrained Set: 3 variables, 2 constraints\n", + " Over-Constrained Set: 0 variables, 0 constraints\n", + "\n", + "------------------------------------------------------------------------------------\n", + "0 Cautions\n", + "\n", + " No cautions found!\n", + "\n", + "------------------------------------------------------------------------------------\n", + "Suggested next steps:\n", + "\n", + " display_underconstrained_set()\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "dt = DiagnosticsToolbox(m)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, the toolbox is reporting one degree of freedom and a structural singularity. In this case, these are both expected as we are trying to solve an optimization problem with one degree of freedom so we can move on to the next step. \n", + "\n", + "### Evaluate the initial point\n", + "\n", + "Before we try to solve our degenerate model, it can often be useful to check the state of the model at the initial step seen by the solver. The cell below creates an instance of Ipopt and sets the iteration limit to 0 so that it terminates at the initial condition." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.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...............: 3.0000000000000000e+00 3.0000000000000000e+00\n", + "Dual infeasibility......: 1.0000000000000000e+00 1.0000000000000000e+00\n", + "Constraint violation....: 2.0000000000000000e+00 2.0000000000000000e+00\n", + "Complementarity.........: 4.0000000499999997e+00 4.0000000499999997e+00\n", + "Overall NLP error.......: 4.0000000499999997e+00 4.0000000499999997e+00\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" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'warning', 'Message': 'Ipopt 3.13.2\\\\x3a Maximum Number of Iterations Exceeded.', 'Termination condition': 'maxIterations', 'Id': 400, 'Error rc': 0, 'Time': 0.004778385162353516}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver = 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", + "solver.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(m2, tee=True)\n", - "\n", - "# Create Degeneracy Hunter object\n", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh2 = DegeneracyHunter(m2, solver=pyo.SolverFactory(\"cbc\"))" + "solver.solve(m, tee=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Identify constraints with residuals greater than 0.1 at the initial point" + "### Identify constraints with large residuals at the initial point\n", + "\n", + "With the solution at the initial point, we can then use the Diagnostics Toolbox to check the residuals of all constraints at the initial point." ] }, { "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "dh2.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 7, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following constraint(s) have large residuals (>1.0E-05):\n", + "\n", + " con2: 2.00000E+00\n", + " con5: 2.00000E+00\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "### Solve the optimization problem and extract the solution" + "# Check for large residuals\n", + "dt.display_constraints_with_large_residuals()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's solve the optimization problem." + "We can also check for variables which are close to (or outside of) their bounds." ] }, { "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2, tee=True)\n", - "\n", - "for i in m2.I:\n", - " print(\"x[\", i, \"]=\", m2.x[i]())" - ] - }, - { - "cell_type": "markdown", + "execution_count": 8, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following variable(s) have values close to their bounds (abs=1.0E-04, rel=1.0E-04):\n", + "\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of linesearch evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution." + "dt.display_variables_near_bounds()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Check the rank of the Jacobian of the equality constraints" + "### Solve the optimization problem and extract the solution\n", + "\n", + "There appear to be no significant issues at the initial point, so let us move on to solving the optimization problem." ] }, { "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "n_deficient = dh2.check_rank_equality_constraints()" - ] - }, - { - "cell_type": "markdown", + "execution_count": 9, "metadata": {}, - "source": [ - "A singular value near 0 indicates the Jacobian of the equality constraints is rank deficient. For each near-zero singular value, there is likely one degenerate constraint." + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 2.17e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01h 1\n", + " 2 1.0184419e+00 1.84e-02 1.26e-02 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.82e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 4.99e-04 - 1.00e+00 5.82e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.33e+01 -2.5 1.88e-04 - 1.00e+00 2.94e-01f 2\n", + " 8 1.0000153e+00 1.53e-05 2.19e+01 -2.5 6.85e-05 - 1.00e+00 8.08e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.78e+02 -2.5 3.47e-05 - 1.00e+00 2.45e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000014e+00 1.44e-06 2.83e-08 -2.5 7.79e-06 - 1.00e+00 1.00e+00h 1\n", + " 11 1.0000000e+00 1.39e-09 1.84e-11 -5.7 2.56e-06 - 1.00e+00 1.00e+00h 1\n", + " 12 1.0000000e+00 1.39e-09 1.24e+02 -8.6 1.16e-08 - 1.00e+00 9.54e-07h 21\n", + " 13 9.9999999e-01 9.82e-09 1.14e-13 -8.6 1.33e-08 - 1.00e+00 1.00e+00s 22\n", + " 14 9.9999999e-01 9.96e-09 6.46e-14 -8.6 2.29e-10 - 1.00e+00 1.00e+00s 22\n", + " 15 9.9999999e-01 9.96e-09 1.82e+02 -9.0 4.00e-11 - 1.00e+00 9.54e-07h 21\n", + " 16 9.9999999e-01 1.00e-08 1.52e-13 -9.0 5.89e-11 - 1.00e+00 1.00e+00s 22\n", + "\n", + "Number of Iterations....: 16\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999000238915e-01 9.9999999000238915e-01\n", + "Dual infeasibility......: 1.5188693260016487e-13 1.5188693260016487e-13\n", + "Constraint violation....: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "Complementarity.........: 9.2217032157601989e-10 9.2217032157601989e-10\n", + "Overall NLP error.......: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 111\n", + "Number of objective gradient evaluations = 17\n", + "Number of equality constraint evaluations = 111\n", + "Number of inequality constraint evaluations = 111\n", + "Number of equality constraint Jacobian evaluations = 17\n", + "Number of inequality constraint Jacobian evaluations = 17\n", + "Number of Lagrangian Hessian evaluations = 16\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: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.04887509346008301}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver.options[\"max_iter\"] = 50\n", + "solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of line search evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution.\n", + "\n", + "The cell below shows the values of the variables at the solution." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 10, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000099975996 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 0.0 : 5 : False : False : Reals\n" + ] + } + ], "source": [ - "### Identify candidate degenerate constraints" + "m.x.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Degeneracy Hunter first identifies candidate degenerate equations." + "### Check the rank of the Jacobian of the equality constraints\n", + "\n", + "One way to check for the presence of degenerate equations is to calculate the rank of the Jacobian matrix of the equality constraints. A singular value near 0 indicates the Jacobian is rank deficient, and for each near-zero singular value there is likely one degenerate constraint.\n", + "\n", + "The Diagnostics Toolbox contains a Singular Value Decomposition (SVD) toolbox which can be used to determine the number of small singular values in the Jacobian for a model as shown below." ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 11, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "\n", + "Number of Singular Values less than 1.0E-06 is 1\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "ds2 = dh2.find_candidate_equations(verbose=True, tee=True)" + "svd = dt.prepare_svd_toolbox()\n", + "svd.display_rank_of_equality_constraints()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Find irreducible degenerate sets (IDS)" + "As can be seen above, there is one singular value less than 1e-6, suggesting we have one degenerate constraint." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Next, Degeneracy Hunter enumerates through the candidate equations. For each candidate equation, Degenerate Hunter solves a MILP to compute the irreducible degenerate set that must contain the candidate equation." + "### Finding Irreducible Degenerate Sets (IDS)\n", + "\n", + "Next, we can use Degeneracy Hunter to identify irreducible degenerate sets; that is, the smallest set of constraints that contain a redundant constraint.\n", + "\n", + "Degeneracy Hunter first identifies candidate degenerate equations by solving a mixed integer linear program (MILP). It then iterates through the candidate equations and solves a second MILP for each to compute the irreducible degenerate set that must contain the candidate equation.\n", + "\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver, users can specify their solver of choice via a configuration argument. The default option is SCIP (https://scipopt.org/) which is an open-source solver.\n", + "
\n", + "\n", + "The cell below shows how to create an instance of Degeneracy Hunter and generate a report of the irreducible degenerate sets." ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 12, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Candidate Equations\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving Candidates MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Irreducible Degenerate Sets\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model to compute irreducible degenerate set.\n", + "Solving MILP 1 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con2.\n", + "Solving MILP 2 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con5.\n", + "====================================================================================\n", + "Irreducible Degenerate Sets\n", + "\n", + " Irreducible Degenerate Set 0\n", + " nu Constraint Name\n", + " 1.0 con2\n", + " -1.0 con5\n", + "\n", + " Irreducible Degenerate Set 1\n", + " nu Constraint Name\n", + " -1.0 con2\n", + " 1.0 con5\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "ids = dh2.find_irreducible_degenerate_sets(verbose=True)" + "dh = dt.prepare_degeneracy_hunter()\n", + "dh.report_irreducible_degenerate_sets()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Reformulate Example 2" + "As can be seen above, Degeneracy Hunter identified two constraints as potentially redundant (the duplicate constraints ``con2`` and ``con5``). Degeneracy Hunter then reports an Irreducible Degenerate Set for each of these constraints, which in the case as identical. More complex models may have partially overlapping IDS's.\n", + "\n", + "These results show us that ``con2`` and ``con5`` are mutually redundant; i.e., we can eliminate one of these with no effect on the model solution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's reformulate the model by skipping/removing the redundant equality constraint:\n", + "### Reformulate Model to Remove Redundant Constraint\n", "\n", - "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", - "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", - "& x_1 + x_2 + x_3 = 1 \\\\\n", - "& x_2 - 2 x_3 \\leq 1 \\\\\n", - "& x_1 + x_3 \\geq 1\n", - "\\end{align*} $$" + "Now let's reformulate the model by removing the redundant equality constraint ``con5`` and solve the reformulated model." ] }, { "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [], - "source": [ - "m2b = example2(with_degenerate_constraint=False)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 13, "metadata": {}, - "source": [ - "### Solve the reformulated model" + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 1\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 6.30e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 1.34e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01f 1\n", + " 2 1.0184419e+00 1.84e-02 9.15e-03 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.88e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 5.00e-04 - 1.00e+00 5.81e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.34e+01 -2.5 1.88e-04 - 1.00e+00 2.93e-01f 2\n", + " 8 1.0000154e+00 1.54e-05 2.26e+01 -2.5 6.89e-05 - 1.00e+00 8.04e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.98e+02 -2.5 3.64e-05 - 1.00e+00 2.33e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000020e+00 2.04e-06 1.33e+02 -2.5 9.72e-06 - 1.00e+00 8.28e-01h 1\n", + " 11 1.0000016e+00 1.61e-06 2.26e+03 -2.5 4.79e-06 - 1.00e+00 2.12e-01f 2\n", + " 12 1.0000002e+00 2.46e-07 8.72e+02 -2.5 1.20e-06 - 1.00e+00 8.47e-01h 1\n", + " 13 1.0000002e+00 1.95e-07 1.71e+04 -2.5 7.02e-07 - 1.00e+00 2.09e-01f 2\n", + " 14 1.0000000e+00 1.54e-08 3.23e+03 -2.5 1.50e-07 - 1.00e+00 9.21e-01h 1\n", + " 15 1.0000000e+00 1.15e-08 9.99e+04 -2.5 6.89e-08 - 1.00e+00 2.54e-01f 2\n", + " 16 1.0000000e+00 2.22e-16 2.83e-08 -2.5 8.21e-09 - 1.00e+00 1.00e+00h 1\n", + " 17 1.0000000e+00 2.22e-16 4.14e-11 -8.6 8.25e-16 - 1.00e+00 1.00e+00 0\n", + "\n", + "Number of Iterations....: 17\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999999999978e-01 9.9999999999999978e-01\n", + "Dual infeasibility......: 4.1425156707686206e-11 4.1425156707686206e-11\n", + "Constraint violation....: 2.2204460492503131e-16 2.2204460492503131e-16\n", + "Complementarity.........: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "Overall NLP error.......: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 23\n", + "Number of objective gradient evaluations = 18\n", + "Number of equality constraint evaluations = 23\n", + "Number of inequality constraint evaluations = 23\n", + "Number of equality constraint Jacobian evaluations = 18\n", + "Number of inequality constraint Jacobian evaluations = 18\n", + "Number of Lagrangian Hessian evaluations = 17\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.002\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 4, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.006506204605102539}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m2 = build_example(include_redundant_constraint=False)\n", + "\n", + "solver.solve(m2, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let us check to see if the optimal solution has changed." ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000005924994 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 4.55844318120227e-11 : 5 : False : False : Reals\n" + ] + } + ], "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2b, tee=True)\n", - "\n", - "for i in m2b.I:\n", - " print(\"x[\", i, \"]=\", m.x[i]())" + "m2.x.display()" ] }, { @@ -650,13 +828,23 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, + "execution_count": 15, + "metadata": { + "tags": [ + "testing" + ] + }, "outputs": [], - "source": [] + "source": [ + "import pytest\n", + "assert pyo.value(m2.x[1]) == pytest.approx(1, rel=1e-6)\n", + "assert pyo.value(m2.x[2]) == pytest.approx(0, abs=1e-6)\n", + "assert pyo.value(m2.x[3]) == pytest.approx(0, abs=1e-6)" + ] } ], "metadata": { + "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", @@ -672,7 +860,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.6" + "version": "3.11.5" } }, "nbformat": 4, diff --git a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_usr.ipynb b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_usr.ipynb index 1b4e368e..ac2388b6 100644 --- a/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_usr.ipynb +++ b/idaes_examples/notebooks/docs/diagnostics/degeneracy_hunter_usr.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": { "tags": [ "header", @@ -34,324 +34,80 @@ "Maintainer: Andrew Lee \n", "Updated: 2023-06-01 \n", "\n", - "This notebook shows how to use the following Degeneracy Hunter features using two motivating examples:\n", + "> Alexander W. Dowling and Lorenz T. Biegler (2015). Degeneracy Hunter: An Algorithm for Determining Irreducible Sets of Degenerate Constraints in Mathematical Programs. 12th International Symposium on Process Systems Engineering and 25th European Symposium on Computer Aided Process Engineering. Ed. by Krist V. Gernaey, Jakob K. Huusom, and Raqul Gani. Computer-Aided Chemical Engineering, Vol. 37, p. 809 \u2013 814. [link to PDF](https://dowlinglab.nd.edu/assets/447116/2015_degeneracy_hunter_an_algorithm_for_determining_irreducible_sets_of_degenerate_constraints_in_mathematical_prog.pdf)\n", + "\n", + "This notebook shows how to use the Degeneracy Hunter tool which is part of the IDAES Diagnostics Toolbox features using two motivating examples:\n", + "\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": "markdown", - "metadata": {}, - "source": [ - "## Setup\n", + "* Demonstrates the performance benefits in Ipopt from removing a single redundant constraint\n", "\n", - "We start by importing Pyomo and Degeneracy Hunter." - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "import pyomo.environ as pyo\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver. This example uses the open-source solver SCIP (https://scipopt.org/), however users may specify another solver of their choice if they have one available.\n", + "
\n", "\n", - "from idaes.core.util.model_diagnostics import DegeneracyHunter" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example 1: Well-Behaved Nonlinear Program" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Consider the following \"well-behaved\" nonlinear optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "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": {}, - "source": [ - "### Define the model in Pyomo" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We start by defining the optimization problem in Pyomo." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "m = pyo.ConcreteModel()\n", + "## What does Degeneracy Mean?\n", "\n", - "m.I = pyo.Set(initialize=[i for i in range(5)])\n", + "In simple terms, a degenerate model contains one or more constraints that do not add any additional information that is not already included in other constraints in the model. The simplest example of this is the case where a constraint in the model is duplicated. For example, consider the model below:\n", "\n", - "m.x = pyo.Var(m.I, bounds=(-10, 10), initialize=1.0)\n", + "$$\\begin{align*}\n", + "& x_1 + x_2 = 1 \\\\\n", + "& x_1 + x_2 = 1 \\\\\n", + "\\end{align*} $$\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", + "Notice the two equality constraints are identical and thus redundant. This means the constraint qualifications that solvers rely on (e.g., [Linear Independence Constraint Qualification or LICQ](https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions#Regularity_conditions_(or_constraint_qualifications))) do not hold which has three important implications:\n", "\n", - "m.obj = pyo.Objective(expr=sum(m.x[i] ** 2 for i in m.I))\n", + "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", + "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", + "3. Theoretical convergence properties of optimization algorithms may not hold\n", "\n", - "m.pprint()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Evaluate the initial point" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": 3, - "metadata": {}, - "outputs": [], - "source": [ - "# Specify Ipopt as the solver\n", - "opt = pyo.SolverFactory(\"ipopt\")\n", + "Put another way, for a deterministic (square) problem we required that the number of equality constraints be equal to the number of free variables. In this case, there appear to be two equality constraints but in reality there is effectively only one as the two constraints are the same. Thus, the problem is not well defined and there exist a family of solutions that can satisfy the constraints given; any combination of values for ``x1`` and ``x2`` that sum to one will satisfy the constraints as written.\n", "\n", - "# Specifying an iteration limit of 0 allows us to inspect the initial point\n", - "opt.options[\"max_iter\"] = 0\n", + "In practice, degeneracies are rarely as straight-forward as a duplicated constraint. A more common occurrence is a set of linearly dependent constraints; that is a system of constraints where one constraint in the set can be formed by a linear combination of other constraints in the set. For example:\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", + "$$\\begin{align*}\n", + "& x_1 + x_2 + x_3 = 1 \\\\\n", + "& 2x_1 + 2x_2 + 2x_3 = 2 \\\\\n", + "\\end{align*} $$\n", "\n", - "# Create Degeneracy Hunter object\n", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh = DegeneracyHunter(m, solver=pyo.SolverFactory(\"cbc\"))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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": {}, - "source": [ - "### Identify the constraint residuals larger than 0.1" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "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", + "Here, the second constraint can be formed by multiplying the first constraint by two. Another common example of linearly dependent constraints is the combination of a set of material balances for all components in a system and an overall material balance. In this case, the overall material balances can be formed by adding all the individual component balances together. Depending on the number of components being modeled, this system of constraints can be quite large, but the end result is the same; the model effectively has one less constraint than it would appear to have.\n", "\n", - "The following line of code will print out all constraints with residuals larger than `0.1`:" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Important: Ipopt does several preprocessing steps when we executed it with zero iterations. When checking the initial point, it is strongly recommended to call Ipopt with zero iterations first. Otherwise, you will not be analyzing the initial point Ipopt starts with." + "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the purpose of Degeneracy Hunter. Let's see it in action." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Identify all variables within 1 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Another common question when developing optimization models is, \"Which variables are within their bounds by a given tolerance?\" Below is the syntax:" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_variable_bounds(tol=1.0)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Solve the optimization problem" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we can solve the optimization problem. We first set the number of iterations to 50 and then resolve with Ipopt." - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m, tee=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, Ipopt has no trouble solving this optimization problem." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check if any constraint residuals are large than 1E-14" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Let's now inspect the new solution to see which (if any) constraints have residuals larger than 10$^{-14}$." - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_residuals(tol=1e-14)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "As expected, all of the constraints are satisfied, even with this fairly tight tolerance." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Identify all variables within 1E-5 of their bounds" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Finally, let's check if any of the variables are near their bounds at the new solution." - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [], - "source": [ - "dh.check_variable_bounds(tol=1e-5)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Great, no variables are near their bounds. If a variable was at its bound, it is important the inspect the model and confirm the bound is physically sensible/what you intended." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Check the rank of the constraint Jacobian at the solution" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The main feature of Degeneracy Hunter is to check if an optimization problem is poorly formulated. Let's see what happens when we check the rank on a carefully formulated optimization problem:" + "## Setup\n", + "\n", + "We start by importing Pyomo and the IDAES Diagnostics Toolbox." ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "dh.check_rank_equality_constraints()" + "import pyomo.environ as pyo\n", + "\n", + "from idaes.core.util.model_diagnostics import DiagnosticsToolbox" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Example 2: Linear Program with Redundant Equality Constraints" + "## Example: Linear Program with Redundant Equality Constraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's apply Degeneracy Hunter to a poorly formulated optimization problem:\n", + "Let's apply these tools to a poorly formulated optimization problem:\n", "\n", "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", @@ -361,256 +117,661 @@ "& x_1 + x_2 + x_3 = 1 \\\\\n", "\\end{align*} $$\n", "\n", + "As you can see, this is similar to our example above with duplicated equality constraints.\n", "\n", - "Notice the two equality constraints are redundant. This means the constraint qualifications (e.g., LICQ) do not hold which has three important implications:\n", - "1. The optimal solution may not be mathematically well-defined (e.g., the dual variables are not unique)\n", - "2. The calculations performed by the optimization solver may become numerically poorly scaled\n", - "3. Theoretical convergence properties of optimization algorithms may not hold\n", - "\n", - "The absolute best defense against this is to detect degenerate equations and reformulate the model to remove them; this is the primary purpose of Degeneracy Hunter. Let's see it in action." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Define the model in Pyomo" + "The cell below shows how to construct this model in Pyomo." ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ - "def example2(with_degenerate_constraint=True):\n", - " \"\"\"Create the Pyomo model for Example 2\n", - "\n", - " Arguments:\n", - " with_degenerate_constraint: Boolean, if True, include the redundant linear constraint\n", - "\n", - " Returns:\n", - " m2: Pyomo model\n", - " \"\"\"\n", - "\n", - " m2 = pyo.ConcreteModel()\n", - "\n", - " m2.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", + "def build_example(include_redundant_constraint=True):\n", + " m = pyo.ConcreteModel()\n", "\n", - " m2.x = pyo.Var(m2.I, bounds=(0, 5), initialize=1.0)\n", + " m.I = pyo.Set(initialize=[i for i in range(1, 4)])\n", "\n", - " m2.con1 = pyo.Constraint(expr=m2.x[1] + m2.x[2] >= 1)\n", - " m2.con2 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", - " m2.con3 = pyo.Constraint(expr=m2.x[2] - 2 * m2.x[3] <= 1)\n", - " m2.con4 = pyo.Constraint(expr=m2.x[1] + m2.x[3] >= 1)\n", + " m.x = pyo.Var(m.I, bounds=(0, 5), initialize=1.0)\n", "\n", - " if with_degenerate_constraint:\n", - " m2.con5 = pyo.Constraint(expr=m2.x[1] + m2.x[2] + m2.x[3] == 1)\n", + " m.con1 = pyo.Constraint(expr=m.x[1] + m.x[2] >= 1)\n", + " m.con2 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", + " m.con3 = pyo.Constraint(expr=m.x[2] - 2 * m.x[3] <= 1)\n", + " m.con4 = pyo.Constraint(expr=m.x[1] + m.x[3] >= 1)\n", + " \n", + " if include_redundant_constraint:\n", + " m.con5 = pyo.Constraint(expr=m.x[1] + m.x[2] + m.x[3] == 1)\n", "\n", - " m2.obj = pyo.Objective(expr=sum(m2.x[i] for i in m2.I))\n", + " m.obj = pyo.Objective(expr=sum(m.x[i] for i in m.I))\n", + " \n", + " return m\n", "\n", - " m2.pprint()\n", - "\n", - " return m2\n", - "\n", - "\n", - "# Create the Pyomo model for Example 2 including the redundant constraint\n", - "m2 = example2()" + "m = build_example()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Evaluate the initial point" + "### Check for Structural Issues\n", + "\n", + "Before we try to solve the model, the first thing we should do is use the Diagnostics Toolbox to check for any structural issues in our model. The cell below creates an instance of the Diagnostics Toolbox and calls the ``report_structural_issues`` method." ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 5, "metadata": {}, - "outputs": [], - "source": [ + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "Model Statistics\n", + "\n", + " Activated Blocks: 1 (Deactivated: 0)\n", + " Free Variables in Activated Constraints: 3 (External: 0)\n", + " Free Variables with only lower bounds: 0\n", + " Free Variables with only upper bounds: 0\n", + " Free Variables with upper and lower bounds: 3\n", + " Fixed Variables in Activated Constraints: 0 (External: 0)\n", + " Activated Equality Constraints: 2 (Deactivated: 0)\n", + " Activated Inequality Constraints: 3 (Deactivated: 0)\n", + " Activated Objectives: 1 (Deactivated: 0)\n", + "\n", + "------------------------------------------------------------------------------------\n", + "2 WARNINGS\n", + "\n", + " WARNING: 1 Degree of Freedom\n", + " WARNING: Structural singularity found\n", + " Under-Constrained Set: 3 variables, 2 constraints\n", + " Over-Constrained Set: 0 variables, 0 constraints\n", + "\n", + "------------------------------------------------------------------------------------\n", + "0 Cautions\n", + "\n", + " No cautions found!\n", + "\n", + "------------------------------------------------------------------------------------\n", + "Suggested next steps:\n", + "\n", + " display_underconstrained_set()\n", + "\n", + "====================================================================================\n" + ] + } + ], + "source": [ + "dt = DiagnosticsToolbox(m)\n", + "dt.report_structural_issues()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you can see, the toolbox is reporting one degree of freedom and a structural singularity. In this case, these are both expected as we are trying to solve an optimization problem with one degree of freedom so we can move on to the next step. \n", + "\n", + "### Evaluate the initial point\n", + "\n", + "Before we try to solve our degenerate model, it can often be useful to check the state of the model at the initial step seen by the solver. The cell below creates an instance of Ipopt and sets the iteration limit to 0 so that it terminates at the initial condition." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.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...............: 3.0000000000000000e+00 3.0000000000000000e+00\n", + "Dual infeasibility......: 1.0000000000000000e+00 1.0000000000000000e+00\n", + "Constraint violation....: 2.0000000000000000e+00 2.0000000000000000e+00\n", + "Complementarity.........: 4.0000000499999997e+00 4.0000000499999997e+00\n", + "Overall NLP error.......: 4.0000000499999997e+00 4.0000000499999997e+00\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" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'warning', 'Message': 'Ipopt 3.13.2\\\\x3a Maximum Number of Iterations Exceeded.', 'Termination condition': 'maxIterations', 'Id': 400, 'Error rc': 0, 'Time': 0.004778385162353516}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver = 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", + "solver.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(m2, tee=True)\n", - "\n", - "# Create Degeneracy Hunter object\n", - "# The Degeneracy Hunter algorithm needs a MILP solver\n", - "# Here we specify CBC, an open source solver\n", - "dh2 = DegeneracyHunter(m2, solver=pyo.SolverFactory(\"cbc\"))" + "solver.solve(m, tee=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Identify constraints with residuals greater than 0.1 at the initial point" + "### Identify constraints with large residuals at the initial point\n", + "\n", + "With the solution at the initial point, we can then use the Diagnostics Toolbox to check the residuals of all constraints at the initial point." ] }, { "cell_type": "code", - "execution_count": 12, - "metadata": {}, - "outputs": [], - "source": [ - "dh2.check_residuals(tol=0.1)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 7, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following constraint(s) have large residuals (>1.0E-05):\n", + "\n", + " con2: 2.00000E+00\n", + " con5: 2.00000E+00\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "### Solve the optimization problem and extract the solution" + "# Check for large residuals\n", + "dt.display_constraints_with_large_residuals()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's solve the optimization problem." + "We can also check for variables which are close to (or outside of) their bounds." ] }, { "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [], - "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2, tee=True)\n", - "\n", - "for i in m2.I:\n", - " print(\"x[\", i, \"]=\", m2.x[i]())" - ] - }, - { - "cell_type": "markdown", + "execution_count": 8, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "The following variable(s) have values close to their bounds (abs=1.0E-04, rel=1.0E-04):\n", + "\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of linesearch evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution." + "dt.display_variables_near_bounds()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Check the rank of the Jacobian of the equality constraints" + "### Solve the optimization problem and extract the solution\n", + "\n", + "There appear to be no significant issues at the initial point, so let us move on to solving the optimization problem." ] }, { "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "n_deficient = dh2.check_rank_equality_constraints()" - ] - }, - { - "cell_type": "markdown", + "execution_count": 9, "metadata": {}, - "source": [ - "A singular value near 0 indicates the Jacobian of the equality constraints is rank deficient. For each near-zero singular value, there is likely one degenerate constraint." + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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...: 6\n", + "Number of nonzeros in inequality constraint Jacobian.: 6\n", + "Number of nonzeros in Lagrangian Hessian.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 2.17e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01h 1\n", + " 2 1.0184419e+00 1.84e-02 1.26e-02 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.82e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 4.99e-04 - 1.00e+00 5.82e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.33e+01 -2.5 1.88e-04 - 1.00e+00 2.94e-01f 2\n", + " 8 1.0000153e+00 1.53e-05 2.19e+01 -2.5 6.85e-05 - 1.00e+00 8.08e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.78e+02 -2.5 3.47e-05 - 1.00e+00 2.45e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000014e+00 1.44e-06 2.83e-08 -2.5 7.79e-06 - 1.00e+00 1.00e+00h 1\n", + " 11 1.0000000e+00 1.39e-09 1.84e-11 -5.7 2.56e-06 - 1.00e+00 1.00e+00h 1\n", + " 12 1.0000000e+00 1.39e-09 1.24e+02 -8.6 1.16e-08 - 1.00e+00 9.54e-07h 21\n", + " 13 9.9999999e-01 9.82e-09 1.14e-13 -8.6 1.33e-08 - 1.00e+00 1.00e+00s 22\n", + " 14 9.9999999e-01 9.96e-09 6.46e-14 -8.6 2.29e-10 - 1.00e+00 1.00e+00s 22\n", + " 15 9.9999999e-01 9.96e-09 1.82e+02 -9.0 4.00e-11 - 1.00e+00 9.54e-07h 21\n", + " 16 9.9999999e-01 1.00e-08 1.52e-13 -9.0 5.89e-11 - 1.00e+00 1.00e+00s 22\n", + "\n", + "Number of Iterations....: 16\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999000238915e-01 9.9999999000238915e-01\n", + "Dual infeasibility......: 1.5188693260016487e-13 1.5188693260016487e-13\n", + "Constraint violation....: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "Complementarity.........: 9.2217032157601989e-10 9.2217032157601989e-10\n", + "Overall NLP error.......: 9.9976108502985994e-09 9.9976108502985994e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 111\n", + "Number of objective gradient evaluations = 17\n", + "Number of equality constraint evaluations = 111\n", + "Number of inequality constraint evaluations = 111\n", + "Number of equality constraint Jacobian evaluations = 17\n", + "Number of inequality constraint Jacobian evaluations = 17\n", + "Number of Lagrangian Hessian evaluations = 16\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: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.04887509346008301}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "solver.options[\"max_iter\"] = 50\n", + "solver.solve(m, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We got lucky here. Ipopt implements several algorithmic and numerical safeguards to handle (mildy) degenerate equations. Nevertheless, notice the last column of the Ipopt output labeled `ls`. This is the number of line search evaluations. For iterations 0 to 11, `ls` is 1, which means Ipopt is taking full steps. For iterations 12 to 16, however, `ls` is greater than 20. This means Ipopt is struggling (a little) to converge to the solution.\n", + "\n", + "The cell below shows the values of the variables at the solution." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 10, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000099975996 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 0.0 : 5 : False : False : Reals\n" + ] + } + ], "source": [ - "### Identify candidate degenerate constraints" + "m.x.display()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Degeneracy Hunter first identifies candidate degenerate equations." + "### Check the rank of the Jacobian of the equality constraints\n", + "\n", + "One way to check for the presence of degenerate equations is to calculate the rank of the Jacobian matrix of the equality constraints. A singular value near 0 indicates the Jacobian is rank deficient, and for each near-zero singular value there is likely one degenerate constraint.\n", + "\n", + "The Diagnostics Toolbox contains a Singular Value Decomposition (SVD) toolbox which can be used to determine the number of small singular values in the Jacobian for a model as shown below." ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 11, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "====================================================================================\n", + "\n", + "Number of Singular Values less than 1.0E-06 is 1\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "ds2 = dh2.find_candidate_equations(verbose=True, tee=True)" + "svd = dt.prepare_svd_toolbox()\n", + "svd.display_rank_of_equality_constraints()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Find irreducible degenerate sets (IDS)" + "As can be seen above, there is one singular value less than 1e-6, suggesting we have one degenerate constraint." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Next, Degeneracy Hunter enumerates through the candidate equations. For each candidate equation, Degenerate Hunter solves a MILP to compute the irreducible degenerate set that must contain the candidate equation." + "### Finding Irreducible Degenerate Sets (IDS)\n", + "\n", + "Next, we can use Degeneracy Hunter to identify irreducible degenerate sets; that is, the smallest set of constraints that contain a redundant constraint.\n", + "\n", + "Degeneracy Hunter first identifies candidate degenerate equations by solving a mixed integer linear program (MILP). It then iterates through the candidate equations and solves a second MILP for each to compute the irreducible degenerate set that must contain the candidate equation.\n", + "\n", + "
\n", + "Note:\n", + "Degeneracy Hunter requires a suitable MILP solver, users can specify their solver of choice via a configuration argument. The default option is SCIP (https://scipopt.org/) which is an open-source solver.\n", + "
\n", + "\n", + "The cell below shows how to create an instance of Degeneracy Hunter and generate a report of the irreducible degenerate sets." ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 12, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Candidate Equations\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving Candidates MILP model.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Searching for Irreducible Degenerate Sets\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Building MILP model to compute irreducible degenerate set.\n", + "Solving MILP 1 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con2.\n", + "Solving MILP 2 of 2.\n", + "2023-11-17 14:33:20 [INFO] idaes.core.util.model_diagnostics: Solving IDS MILP for constraint con5.\n", + "====================================================================================\n", + "Irreducible Degenerate Sets\n", + "\n", + " Irreducible Degenerate Set 0\n", + " nu Constraint Name\n", + " 1.0 con2\n", + " -1.0 con5\n", + "\n", + " Irreducible Degenerate Set 1\n", + " nu Constraint Name\n", + " -1.0 con2\n", + " 1.0 con5\n", + "\n", + "====================================================================================\n" + ] + } + ], "source": [ - "ids = dh2.find_irreducible_degenerate_sets(verbose=True)" + "dh = dt.prepare_degeneracy_hunter()\n", + "dh.report_irreducible_degenerate_sets()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Reformulate Example 2" + "As can be seen above, Degeneracy Hunter identified two constraints as potentially redundant (the duplicate constraints ``con2`` and ``con5``). Degeneracy Hunter then reports an Irreducible Degenerate Set for each of these constraints, which in the case as identical. More complex models may have partially overlapping IDS's.\n", + "\n", + "These results show us that ``con2`` and ``con5`` are mutually redundant; i.e., we can eliminate one of these with no effect on the model solution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's reformulate the model by skipping/removing the redundant equality constraint:\n", + "### Reformulate Model to Remove Redundant Constraint\n", "\n", - "$$\\begin{align*} \\min_{\\mathbf{x}} \\quad & \\sum_{i=\\{1,...,3\\}} x_i \\\\\n", - "\\mathrm{s.t.}~~& x_1 + x_2 \\geq 1 \\\\\n", - "& x_1 + x_2 + x_3 = 1 \\\\\n", - "& x_2 - 2 x_3 \\leq 1 \\\\\n", - "& x_1 + x_3 \\geq 1\n", - "\\end{align*} $$" + "Now let's reformulate the model by removing the redundant equality constraint ``con5`` and solve the reformulated model." ] }, { "cell_type": "code", - "execution_count": 17, - "metadata": {}, - "outputs": [], - "source": [ - "m2b = example2(with_degenerate_constraint=False)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 13, "metadata": {}, - "source": [ - "### Solve the reformulated model" + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt 3.13.2: max_iter=50\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.............: 0\n", + "\n", + "Total number of variables............................: 3\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 3\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 1\n", + "Total number of inequality constraints...............: 3\n", + " inequality constraints with only lower bounds: 2\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 1\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 3.0000000e+00 2.00e+00 6.30e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 1.3934810e+00 3.93e-01 1.34e-01 -1.0 1.23e+00 - 7.90e-01 8.03e-01f 1\n", + " 2 1.0184419e+00 1.84e-02 9.15e-03 -1.7 7.02e-01 - 9.41e-01 9.53e-01h 1\n", + " 3 1.0011246e+00 1.12e-03 3.88e-02 -2.5 1.50e-02 - 1.00e+00 9.39e-01h 1\n", + " 4 1.0006914e+00 6.91e-04 3.12e+00 -2.5 3.17e-03 - 1.00e+00 3.85e-01h 1\n", + " 5 1.0002664e+00 2.66e-04 4.35e+00 -2.5 1.12e-03 - 1.00e+00 6.15e-01h 1\n", + " 6 1.0001115e+00 1.12e-04 1.07e+01 -2.5 5.00e-04 - 1.00e+00 5.81e-01h 1\n", + " 7 1.0000788e+00 7.88e-05 4.34e+01 -2.5 1.88e-04 - 1.00e+00 2.93e-01f 2\n", + " 8 1.0000154e+00 1.54e-05 2.26e+01 -2.5 6.89e-05 - 1.00e+00 8.04e-01h 1\n", + " 9 1.0000118e+00 1.18e-05 2.98e+02 -2.5 3.64e-05 - 1.00e+00 2.33e-01f 2\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 1.0000020e+00 2.04e-06 1.33e+02 -2.5 9.72e-06 - 1.00e+00 8.28e-01h 1\n", + " 11 1.0000016e+00 1.61e-06 2.26e+03 -2.5 4.79e-06 - 1.00e+00 2.12e-01f 2\n", + " 12 1.0000002e+00 2.46e-07 8.72e+02 -2.5 1.20e-06 - 1.00e+00 8.47e-01h 1\n", + " 13 1.0000002e+00 1.95e-07 1.71e+04 -2.5 7.02e-07 - 1.00e+00 2.09e-01f 2\n", + " 14 1.0000000e+00 1.54e-08 3.23e+03 -2.5 1.50e-07 - 1.00e+00 9.21e-01h 1\n", + " 15 1.0000000e+00 1.15e-08 9.99e+04 -2.5 6.89e-08 - 1.00e+00 2.54e-01f 2\n", + " 16 1.0000000e+00 2.22e-16 2.83e-08 -2.5 8.21e-09 - 1.00e+00 1.00e+00h 1\n", + " 17 1.0000000e+00 2.22e-16 4.14e-11 -8.6 8.25e-16 - 1.00e+00 1.00e+00 0\n", + "\n", + "Number of Iterations....: 17\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 9.9999999999999978e-01 9.9999999999999978e-01\n", + "Dual infeasibility......: 4.1425156707686206e-11 4.1425156707686206e-11\n", + "Constraint violation....: 2.2204460492503131e-16 2.2204460492503131e-16\n", + "Complementarity.........: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "Overall NLP error.......: 2.6658012082875325e-09 2.6658012082875325e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 23\n", + "Number of objective gradient evaluations = 18\n", + "Number of equality constraint evaluations = 23\n", + "Number of inequality constraint evaluations = 23\n", + "Number of equality constraint Jacobian evaluations = 18\n", + "Number of inequality constraint Jacobian evaluations = 18\n", + "Number of Lagrangian Hessian evaluations = 17\n", + "Total CPU secs in IPOPT (w/o function evaluations) = 0.002\n", + "Total CPU secs in NLP function evaluations = 0.000\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 4, 'Number of variables': 3, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.006506204605102539}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m2 = build_example(include_redundant_constraint=False)\n", + "\n", + "solver.solve(m2, tee=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let us check to see if the optimal solution has changed." ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x : Size=3, Index=I\n", + " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", + " 1 : 0 : 1.0000000005924994 : 5 : False : False : Reals\n", + " 2 : 0 : 0.0 : 5 : False : False : Reals\n", + " 3 : 0 : 4.55844318120227e-11 : 5 : False : False : Reals\n" + ] + } + ], "source": [ - "opt.options[\"max_iter\"] = 50\n", - "opt.solve(m2b, tee=True)\n", - "\n", - "for i in m2b.I:\n", - " print(\"x[\", i, \"]=\", m.x[i]())" + "m2.x.display()" ] }, { @@ -647,16 +808,10 @@ "\n", "Often degenerate equations have a much worse impact on large-scale problems; for example, degenerate equations can cause Ipopt to require many more iterations or terminate at an infeasible point.\n" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { + "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", @@ -672,7 +827,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.6" + "version": "3.11.5" } }, "nbformat": 4,