{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear Optimization in Julia" ], "id": "d8f72118-5fd1-4494-b31b-db4737bb80b2" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using Plots\n", "using JuMP\n", "using HiGHS" ], "id": "4" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overview\n", "\n", "This tutorial will demonstrate how to solve linear optimization problems\n", "graphically and using the `JuMP` package in Julia. It draws heavily from\n", "[this\n", "tutorial](https://github.com/Power-Systems-Optimization-Course/power-systems-optimization/blob/master/Notebooks/02-Anatomy-of-a-Model.ipynb)\n", "by Jesse D. Jenkins and Michael R. Davidson.\n", "\n", "[JuMP](https://jump.dev/) (“**Ju**lia for **M**athematical\n", "**P**rogramming”) is an open-source Julia package that adds\n", "functionality for formulating and solving a variety of optimization\n", "problems. One advantage of JuMP is that its syntax matches the typical\n", "mathematical formalism used to specify optimization problems. We will\n", "use JuMP in this class for our optimization work.\n", "\n", "> **Read the Documentation!**\n", ">\n", "> Make sure that you take a look at the [`JuMP`\n", "> documentation](https://jump.dev/JuMP.jl/stable/) whenever you have a\n", "> question or want to find out how to do something that we don’t discuss\n", "> in any of our tutorials or lectures (or how to do it better!).\n", "\n", "## Setup\n", "\n", "Here we will outline the basic steps for configuring JuMP, though you\n", "can also refer to the official [Installation\n", "Guide](https://jump.dev/JuMP.jl/stable/installation/#Installation-Guide).\n", "\n", "If `JuMP` is not already in your environment (it will be for any of your\n", "assignments, but may not be if you’re doing something independently),\n", "you will need to install it. You will also need to [select a solver and\n", "install the relevant\n", "package](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers).\n", "Some of these are commercial, while others are open source. Solvers are\n", "also not typically universal, as different types of optimization\n", "problems use different algorithms, so be aware of what problem you’re\n", "trying to solve instead of just blindly copying code from one task to\n", "another.\n", "\n", "For example, for the linear programming example, we will use the [HiGHS\n", "solver](https://highs.dev/) via the\n", "[`HiGHS.jl`](https://github.com/jump-dev/HiGHS.jl) package. As seen on\n", "the solver table, HiGHS is open source (via the [MIT\n", "license](https://choosealicense.com/licenses/mit/)) and can solve linear\n", "programs (LP) and mixed-integer linear programs (MILP), as well as\n", "quadratic programs (which we won’t discuss in this course).\n", "\n", "## Linear Programming Example: How Many Widgets Should A Factory Produce?\n", "\n", "### Defining The Problem\n", "\n", "Suppose we own a factory that can produce two types of widgets:\n", "\n", "- Widget A generates a profit of $p_A = \\$100$ per widget; and\n", "- Widget B generates a profit of $p_B = \\$75$ per widget.\n", "\n", "Let $x$ be the number of units of widget A that we want to produce, and\n", "$y$ the number of units of widget B. Our goal is to *maximize* our total\n", "profit $p_Ax + p_By$. This is the **objective function**. We express\n", "this objective using the equation\n", "\n", "This isn’t a very interesting problem yet! We would simply build as much\n", "of both widgets as we could, because there are no constraints on our\n", "ability to produce. To make this more realistic, let’s suppose that both\n", "widgets are produced using the same raw material $M$, of which we can\n", "only procure 300 units. Then, if:\n", "\n", "- Widget A requires 40 units of $M$ per widget, and\n", "- Widget B requires 20 units of $M$ per widget, we arrive at the\n", " following material **constraint**:\n", "\n", "But we might have another constraint: time! Each widget may take a\n", "different amount of labor to produce. For example, let’s say that\n", "\n", "- Widget A takes 6 hours to produce and\n", "- Widget B takes 12 hours to produce. Further, there are only 80 hours\n", " per work that can be allocated to widget production. This becomes the\n", " time constraint\n", "\n", "Finally, we cannot build a negative number of either type of widget.\n", "This is known as a *non-negativity constraint*, and can be expressed as\n", "\n", "Consolidating equations – gives us the following **constrained\n", "optimization problem**:\n", "\n", "### Visualizing the Problem\n", "\n", "Let’s do some plotting to examine the geometry of our optimization\n", "problem. We can do this using the\n", "[`Plots.jl`](https://docs.juliaplots.org/latest/) package in Julia." ], "id": "d3a83d1c-a06e-4fd9-a2b8-ed20ba7f108e" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "## set up objective function parameters and variables\n", "a = range(0, 8, step=0.25)\n", "b = range(0, 8, step=0.25)\n", "\n", "## define objective function\n", "pa = 100\n", "pb = 75\n", "f(a, b) = pa * a + pb * b\n", "\n", "## start plotting\n", "contour(a,b,(a,b)->f(a,b),nlevels=15, c=:heat, linewidth=10, colorbar = false, contour_labels = true) # objective function contours\n", "title!(\"Factory Optimization Problem\") # add title\n", "xaxis!(\"x=Widget A\", lims=(0, maximum(a))) # add x-axis title and limits\n", "yaxis!(\"y=Widget B\", lims=(0, maximum(b))) # add y-axis title and limits\n", "xticks!(0:maximum(a)) # set x-axis ticks\n", "yticks!(0:maximum(b)) # set y-axis ticks\n", "areaplot!(a[a.<=11], (300 .- 40*a)./20, legend=false, opacity=0.3) # plot materials constraint feasible region\n", "areaplot!(a[a.<=8], (80 .- 6*a)./12, legend=false, opacity=0.3) # plot time constraint feasible region" ], "id": "cell-fig-widget-space" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see exactly where the solution will be in\n", "Figure 1, at the\n", "intersection of the feasible regions imposed by the two constraints!\n", "\n", "> **Objective Function Gradient and Solution Uniqueness**\n", ">\n", "> What would happen if one of the constraints were parallel to the level\n", "> sets of the objective function?\n", "\n", "Let’s now use JuMP to identify the location of this point (though we\n", "could also solve for it using linear algebra).\n", "\n", "## Solving This Problem Using JuMP\n", "\n", "### Setting Up the Model and Solver\n", "\n", "To solve our problem, first we need to **define the model**. The *model\n", "object* has lots of attributes, including the variables, constraints,\n", "solver options, etc. We create a new model using the `Model()` function.\n", "Since we are using the `HiGHS` solver, we need to tell JuMP to use the\n", "`HiGHS.Optimizer` solver function." ], "id": "bb71cb3a-7c65-47dc-b5c1-7540ef873a70" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "factory_model = Model(HiGHS.Optimizer)" ], "id": "8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are a bunch of attributes and options that we could set, but we\n", "won’t in this example. If needed, look at the `HiGHS.jl`\n", "[documentation](https://github.com/jump-dev/HiGHS.jl).\n", "\n", "#### Define Variables\n", "\n", "Decision variables ($x$ and $y$ in this case) in JuMP are defined using\n", "the `@variable` macro. The first argument passed to `@variable()` is the\n", "model object, in this case, `factory_model`, and the second argument are\n", "bounds on that variable, created using `>=` and `<=`. JuMP will\n", "interpret the bound specification to obtain the variable name. In this\n", "case, our only bounds directly on the variables are the non-negativity\n", "constraints." ], "id": "e669fa1e-7bf9-4d34-92ed-82867b8d3aed" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "@variable(factory_model, x >= 0)\n", "@variable(factory_model, y >= 0)" ], "id": "10" }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we had a free (or unbounded) variable $z$, we could declare that\n", "variable using `@variable(model, z)`. JuMP also requires unique names\n", "for each variable, or it will throw an error. This is one place where\n", "it’s nice that Julia lets us use sub- and superscripts in variable\n", "names!\n", "\n", "If we did want to modify the bounds after defining the variable, we\n", "could do so using the `set_lower_bound()` and `set_upper_bound()`\n", "functions, or we could remove them using `delete_lower_bound()` and\n", "`delete_upper_bound()`.\n", "\n", "Finally, if we want to see all of the variables associated with a model,\n", "we can use the `all_variables` function to obtain an array." ], "id": "444f5ad3-62a8-449a-ba20-37a8602a5ffe" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "2-element Vector{VariableRef}:\n", " x\n", " y" ] } } ], "source": [ "all_variables(factory_model)" ], "id": "12" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define Constraints\n", "\n", "When defining variables, we were able to declare constraints on their\n", "values by specifying upper and lower bounds. However, we also have other\n", "constraints, which involve multiple decision variables. These are\n", "specified using the `@constraint` macro. We also can pass names for each\n", "constraint[1]. We will use `time` for the time constraint and\n", "`materials` for the materials constraint. These names must be unique.\n", "\n", "[1] Constraint names aren’t required, but are useful if you want to\n", "modify constraints or [get dual variables/shadow\n", "prices](http://localhost:4200/tutorials/julia-jump.html#dual-solutions)." ], "id": "3021d8b0-2b9d-4013-95da-b84b5f697e0e" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "@constraint(factory_model, time, 6x + 12y <= 80) # specify the time constraint\n", "@constraint(factory_model, materials, 40x + 20y <= 300) # materials constraint" ], "id": "14" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define Objective Function\n", "\n", "So far, we’ve defined the feasible region of the decision-variable\n", "domain by setting the constraints. But we need to specify our objective\n", "function to know what we are trying to minimize or maximize over this\n", "region. We define the objective function using the `@objective` macro.\n", "In addition to specifying the model objective and the function, we need\n", "to tell JuMP whether we want to minimize or maximize." ], "id": "7178f772-f022-4c4c-9913-152c3beec529" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "@objective(factory_model, Max, 100x + 75y)" ], "id": "16" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Looking At The Full Model\n", "\n", "Now, let’s look at the model specification. `print()` will print out a\n", "formatted version of the model; in a notebook (or on this page), that\n", "will be marked up with LaTeX, in a REPL terminal, it will not be." ], "id": "03f153a1-943a-41ac-9c3b-346d7fa8d328" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Max 100 x + 75 y\n", "Subject to\n", " time : 6 x + 12 y ≤ 80\n", " materials : 40 x + 20 y ≤ 300\n", " x ≥ 0\n", " y ≥ 0" ] } ], "source": [ "print(factory_model)" ], "id": "18" }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want a LaTeX-marked up version in the REPL, use\n", "`latex_formulation()`.\n", "\n", "We won’t go into detail here, but there are other ways to define the\n", "model, which are detailed in the [JuMP\n", "documentation](https://jump.dev/JuMP.jl/stable/). For example, we can\n", "specify multiple variables using\n", "[`@variables`](https://jump.dev/JuMP.jl/stable/manual/variables/#variables).\n", "Similarly, we can use\n", "[`@constraints`](https://jump.dev/JuMP.jl/stable/manual/constraints/#The-@constraints-macro)\n", "to define multiple constraints at once. Or we can use\n", "[loops](https://jump.dev/JuMP.jl/stable/tutorials/Getting%20started/variables_constraints_objective/#Constraints-in-a-loop)\n", "to define multiple constraints or constraints involving many variables.\n", "We can also specify the model in [vectorized\n", "syntax](https://jump.dev/JuMP.jl/stable/tutorials/Getting%20started/variables_constraints_objective/#Vectorized-syntax),\n", "which is similar to how linear programs are specified in MATLAB.\n", "\n", "### Solve the Model\n", "\n", "Now it’s time to solve the model and find the optimal values\n", "$(x^*, y^*)$. Since we specified the solver when we initialized\n", "`factory_model`, all we have to do is call the `optimize!` function." ], "id": "7baccb52-f485-49c0-955e-7508d7a017e1" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Running HiGHS 1.11.0 (git hash: 364c83a51e): Copyright (c) 2025 HiGHS under MIT licence terms\n", "LP has 2 rows; 2 cols; 4 nonzeros\n", "Coefficient ranges:\n", " Matrix [6e+00, 4e+01]\n", " Cost [8e+01, 1e+02]\n", " Bound [0e+00, 0e+00]\n", " RHS [8e+01, 3e+02]\n", "Presolving model\n", "2 rows, 2 cols, 4 nonzeros 0s\n", "2 rows, 2 cols, 4 nonzeros 0s\n", "Presolve : Reductions: rows 2(-0); columns 2(-0); elements 4(-0) - Not reduced\n", "Problem not reduced by presolve: solving the LP\n", "Using EKK dual simplex solver - serial\n", " Iteration Objective Infeasibilities num(sum)\n", " 0 -1.0937489944e+01 Ph1: 2(4.125); Du: 2(10.9375) 0s\n", " 2 8.4722222222e+02 Pr: 0(0) 0s\n", "Model status : Optimal\n", "Simplex iterations: 2\n", "Objective value : 8.4722222222e+02\n", "P-D objective error : 6.7054298414e-17\n", "HiGHS run time : 0.00" ] } ], "source": [ "optimize!(factory_model)" ], "id": "20" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Query the Solution\n", "\n", "To find the optimal values of our decision variables, we need to query\n", "the values of the variables using `value.()`. We use `value.()` (the\n", "vectorized version of `value()`) because JuMP stores decision variables\n", "differently depending on their number and how they were defined. Uses\n", "the dot-syntax here works with any model specification, while the plain\n", "`value()` will not work if a queried decision variable is stored as a\n", "vector." ], "id": "ab900983-a9ab-45d1-8457-f736a38a1a65" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "5.555555555555554" ] } } ], "source": [ "value.(x)" ], "id": "22" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "3.8888888888888897" ] } } ], "source": [ "value.(y)" ], "id": "24" }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we can see that our optimal inputs are $$(x^*, y^*) = (5.56, 3.89)$$\n", "(and we’ll pretend that we can manufacture and sell parts of widgets).\n", "\n", "### Visualize the Solution\n", "\n", "Let’s take our previous plot and add the solution point to make sure\n", "that we got the solution we expected." ], "id": "a03a4326-1f9c-499b-ac2f-cb92b1c1cb31" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "contour(a,b,(a,b)->f(a,b),nlevels=15, c=:heat, linewidth=10, colorbar = false, contour_labels = true) # objective function contours\n", "title!(\"Factory Optimization Problem\") # add title\n", "xaxis!(\"x=Widget A\", lims=(0, maximum(a))) # add x-axis title and limits\n", "yaxis!(\"y=Widget B\", lims=(0, maximum(b))) # add y-axis title and limits\n", "xticks!(0:maximum(a)) # set x-axis ticks\n", "yticks!(0:maximum(b)) # set y-axis ticks\n", "areaplot!(a, (300 .- 40*a)./20, legend=false, opacity=0.3) # plot materials constraint feasible region\n", "areaplot!(a, (80 .- 6*a)./12, legend=false, opacity=0.3) # plot time constraint feasible region\n", "\n", "## now we plot the solution that we obtained\n", "scatter!([value.(x)],[value.(y)], markercolor=\"blue\", markersize=5)" ], "id": "cell-fig-widget-solution" }, { "cell_type": "markdown", "metadata": {}, "source": [ "As shown in\n", "Figure 2, the\n", "optimal solution $$(x^*, y^*)$$ is exactly where we deduced it would be\n", "geometrically.\n", "\n", "### Other Stuff We Can Do\n", "\n", "We can also use `value.()` to evaluate our constraints without manually\n", "using the equations." ], "id": "3aae9cf8-6693-49c0-b78a-4e87bccba888" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "80.0" ] } } ], "source": [ "value.(time)" ], "id": "28" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "300.0" ] } } ], "source": [ "value.(materials)" ], "id": "30" }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we also want the optimal objective value? We can obtain this\n", "using `objective_value()`." ], "id": "147efa1a-e1d6-49a0-a5ba-d74fca0376c7" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "847.2222222222221" ] } } ], "source": [ "objective_value(factory_model)" ], "id": "32" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could also define other expressions via the [`@expression`\n", "macro](https://jump.dev/JuMP.jl/stable/expressions/#JuMP.@expression) as\n", "functions of the decision variables and evaluate those. For example,\n", "let’s say that we wanted to know the total number of widgets we’d\n", "produce under our optimal allocation of resources." ], "id": "22c1b4b8-d6d6-48b2-80a7-ad9b32e215c6" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "9.444444444444443" ] } } ], "source": [ "@expression(factory_model, total_widgets, x+y)\n", "value.(total_widgets)" ], "id": "34" }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Dual Solutions\n", "\n", "We can identify if our model has a dual solution by calling\n", "`has_duals()`." ], "id": "5d424e7d-d671-4a03-877e-d22b22d7bbd8" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "true" ] } } ], "source": [ "has_duals(factory_model)" ], "id": "36" }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we want to know the dual solution associated with a constraint, we\n", "use the `shadow_price()` function.\n", "\n", "> **Naming Constraints**\n", ">\n", "> Modifying constraints and querying for dual solutions is the reason\n", "> why it’s important to give constraints individual names." ], "id": "70e230bb-a71e-4ac8-a349-58b82eb66936" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "2.7777777777777772" ] } } ], "source": [ "shadow_price(time)" ], "id": "38" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "2.0833333333333335" ] } } ], "source": [ "shadow_price(materials)" ], "id": "40" }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the binding constraint was a variable bound, we could also query that\n", "shadow price by calling `reduced_cost()` on the variable." ], "id": "9cf53ed7-75d5-477d-952d-db98e1876a16" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "-0.0" ] } } ], "source": [ "reduced_cost(x)" ], "id": "42" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "-0.0" ] } } ], "source": [ "reduced_cost(y)" ], "id": "44" }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, the relevant shadow prices are zero because the optimum is\n", "in the interior of the domain. If we had added a strong enough upper\n", "bound on the value(s) of one or both of our decision variables (say,\n", "$x \\leq 4$), then this would be non-zero." ], "id": "0acdc782-c1fa-4c07-9b8a-161d9f452c6d" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernel_info": { "name": "julia" }, "kernelspec": { "name": "julia", "display_name": "Julia", "language": "julia" }, "language_info": { "name": "julia", "codemirror_mode": "julia", "version": "1.11.5" } } }