Skip to content

Instantly share code, notes, and snippets.

@gidden
Created January 22, 2020 08:00
Show Gist options
  • Save gidden/39eb170ca29bb14acee6ea092b9255c3 to your computer and use it in GitHub Desktop.
Save gidden/39eb170ca29bb14acee6ea092b9255c3 to your computer and use it in GitHub Desktop.
Pyomo Differentiate Example
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from sympy import symbols, diff\n",
"from sympy.utilities.lambdify import lambdify\n",
"\n",
"from pyomo.environ import ConcreteModel, PercentFraction, Var, Objective, Constraint\n",
"from pyomo.opt import SolverFactory"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example Problem\n",
"\n",
"\\begin{align}\n",
"max: 3 x_1 + x_2 \\\\\n",
"s.t. y = x_1^2 + x_2^2 \\\\\n",
"y \\leq 0.5 \\\\\n",
"\\frac{\\partial y}{\\partial x_1} \\leq \\frac{1}{2} \\\\\n",
"x_1 \\in [0, 1] \\\\\n",
"x_2 \\in [0, 1] \\\\\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def solve(model):\n",
" solver = 'ipopt'\n",
" solver_io = 'nl'\n",
" opt = SolverFactory(solver, solver_io=solver_io)\n",
" opt.solve(model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sympy Implementation"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1.4114378434226194,\n",
" 0.2500000040426455,\n",
" 0.2500000040426455,\n",
" 0.5000000066862389)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x_1, x_2 = symbols('x_1 x_2')\n",
"y = x_1 ** 2 + x_2 **2\n",
"\n",
"model = ConcreteModel()\n",
"\n",
"model.x = Var([0, 1], domain=PercentFraction)\n",
"model.y = Var([0])\n",
"\n",
"model.OBJ = Objective(expr=3 * model.x[0] + model.x[1], sense=-1)\n",
"model.Constraint1 = Constraint(expr=model.y[0] == model.x[0] ** 2 + model.x[1] ** 2)\n",
"model.Constraint2 = Constraint(expr=model.y[0] <= 0.5)\n",
"model.Constraint3 = Constraint(expr=lambdify(x_1, diff(y, x_1))(model.x[0]) <= 0.5)\n",
"\n",
"solve(model)\n",
"model.OBJ(), model.x[0].value, model.x[0].value, model.y[0].value"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Differentiate Implementation (incorrect)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"from pyomo.core.expr import differentiate"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ERROR: Constructing component 'Constraint3' from data=None failed: ValueError:\n",
" Invalid constraint expression. The constraint expression resolved to a\n",
" trivial Boolean (True) instead of a Pyomo object. Please modify your rule\n",
" to return Constraint.Feasible instead of True.\n",
"\n",
" Error thrown for Constraint 'Constraint3'\n"
]
},
{
"ename": "ValueError",
"evalue": "Invalid constraint expression. The constraint expression resolved to a trivial Boolean (True) instead of a Pyomo object. Please modify your rule to return Constraint.Feasible instead of True.\n\nError thrown for Constraint 'Constraint3'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-5-dfd202bbbaf6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mConstraint1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mConstraint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m**\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m**\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mConstraint2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mConstraint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mConstraint3\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mConstraint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdifferentiate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mConstraint1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbody\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwrt\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 11\u001b[0m \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/work/miniconda3/envs/py37/lib/python3.7/site-packages/pyomo/core/base/block.py\u001b[0m in \u001b[0;36m__setattr__\u001b[0;34m(self, name, val)\u001b[0m\n\u001b[1;32m 568\u001b[0m \u001b[0;31m# Pyomo components are added with the add_component method.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 569\u001b[0m \u001b[0;31m#\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 570\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_component\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mval\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 571\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 572\u001b[0m \u001b[0;31m#\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/work/miniconda3/envs/py37/lib/python3.7/site-packages/pyomo/core/base/block.py\u001b[0m in \u001b[0;36madd_component\u001b[0;34m(self, name, val)\u001b[0m\n\u001b[1;32m 1008\u001b[0m _blockName, str(data))\n\u001b[1;32m 1009\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1010\u001b[0;31m \u001b[0mval\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconstruct\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1011\u001b[0m \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1012\u001b[0m \u001b[0merr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexc_info\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/work/miniconda3/envs/py37/lib/python3.7/site-packages/pyomo/core/base/constraint.py\u001b[0m in \u001b[0;36mconstruct\u001b[0;34m(self, data)\u001b[0m\n\u001b[1;32m 763\u001b[0m err))\n\u001b[1;32m 764\u001b[0m \u001b[0;32mraise\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 765\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_setitem_when_not_present\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtmp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 766\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 767\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/work/miniconda3/envs/py37/lib/python3.7/site-packages/pyomo/core/base/constraint.py\u001b[0m in \u001b[0;36m_setitem_when_not_present\u001b[0;34m(self, index, value)\u001b[0m\n\u001b[1;32m 712\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 713\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_setitem_when_not_present\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 714\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_check_skip_add\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 715\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 716\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/work/miniconda3/envs/py37/lib/python3.7/site-packages/pyomo/core/base/constraint.py\u001b[0m in \u001b[0;36m_check_skip_add\u001b[0;34m(self, index, expr)\u001b[0m\n\u001b[1;32m 889\u001b[0m \u001b[0mexpr\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;34m\"Feasible\"\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m\"Infeasible\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 890\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 891\u001b[0;31m _get_indexed_component_data_name(self,index)))\n\u001b[0m\u001b[1;32m 892\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 893\u001b[0m \u001b[0;31m#\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: Invalid constraint expression. The constraint expression resolved to a trivial Boolean (True) instead of a Pyomo object. Please modify your rule to return Constraint.Feasible instead of True.\n\nError thrown for Constraint 'Constraint3'"
]
}
],
"source": [
"model = ConcreteModel()\n",
"\n",
"model.x = Var([0, 1], domain=PercentFraction, initialize=0.3)\n",
"model.y = Var([0], domain=PercentFraction, initialize=0.4)\n",
"\n",
"model.OBJ = Objective(expr=3 * model.x[0] + model.x[1], sense=-1)\n",
"model.Constraint1 = Constraint(expr=model.y[0] == model.x[0] ** 2 + model.x[1] ** 2)\n",
"model.Constraint2 = Constraint(expr=model.y[0] <= 0.5)\n",
"model.Constraint3 = Constraint(expr=differentiate(model.Constraint1.body, wrt=model.x[0]) <= 0.5)\n",
"\n",
"solve(model)\n",
"model.OBJ(), model.x[0].value, model.x[1].value, model.y[0].value"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment