From 8b86774f1bdfc27897e59cf86d5d0ecfd8572d61 Mon Sep 17 00:00:00 2001 From: Dimitri Alston <123396563+DimitriAlston@users.noreply.github.com> Date: Mon, 22 Jun 2026 16:02:31 -0400 Subject: [PATCH 1/3] Add InfiniteOpt extension --- Project.toml | 9 +- README.md | 5 + examples/infinite_model.ipynb | 422 ++++++++++++++++++++++++++++++++++ examples/infinite_model.jl | 89 +++++++ ext/EOIInfiniteOptExt.jl | 65 ++++++ test/runtests.jl | 60 +++++ 6 files changed, 649 insertions(+), 1 deletion(-) create mode 100644 examples/infinite_model.ipynb create mode 100644 examples/infinite_model.jl create mode 100644 ext/EOIInfiniteOptExt.jl diff --git a/Project.toml b/Project.toml index 039a3a1..fd4a3e0 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,12 @@ SciCompDSL = "91a8cdf1-4ca6-467b-a780-87fda3fff15e" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +[weakdeps] +InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" + +[extensions] +EOIInfiniteOptExt = "InfiniteOpt" + [compat] DocStringExtensions = "0.9" JuMP = "1" @@ -25,8 +31,9 @@ julia = "1.10" [extras] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["CSV", "DataFrames", "Ipopt", "Test"] +test = ["CSV", "DataFrames", "InfiniteOpt", "Ipopt", "Test"] diff --git a/README.md b/README.md index 010fa7b..7977fd9 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,10 @@ An optimal reactor-separator-recycle process design problem originally presented A nonlinear kinetic parameter estimation problem originally described by [[5](#references)] is used to demonstrate the use of `register_odesystem` to formulate and solve an ODE system using [`Ipopt`](https://github.com/coin-or/ipopt) [[6](#references)]. +### [InfiniteOpt ODE System](https://github.com/PSORLab/EOptInterface.jl/blob/main/examples/infinite_model.jl) + +The same ODE example as above is formulated and solved using [`InfiniteOpt.jl`](https://github.com/infiniteopt/InfiniteOpt.jl) [[7](#references)]. + ## References 1. Ma, Y., Gowda, S., Anantharaman, R., Laughman, C., Shah, V., and Rackauckas, C. ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling. (2022). DOI: [10.48550/arXiv.2103.05244](https://doi.org/10.48550/arXiv.2103.05244) 2. Lubin, M., Dowson, O., Dias Garcia, J., Huchette, J., Legat, B., and Vielma, J.P. JuMP 1.0: recent improvements to a modeling language for mathematical optimization. *Mathematical Programming Computation.* 15, 581-589 (2023). DOI: [10.1007/s12532-023-00239-3](https://doi.org/10.1007/s12532-023-00239-3) @@ -29,3 +33,4 @@ A nonlinear kinetic parameter estimation problem originally described by [[5](#r 4. Wilhelm, M. E. and Stuber, M.D. EAGO.jl: easy advanced global optimization in Julia. *Optimization Methods & Software.* 37(2), 425-450 (2022). DOI: [10.1080/10556788.2020.1786566](https://doi.org/10.1080/10556788.2020.1786566) 5. Taylor, J.W. Direct measurement and analysis of cyclohexadienyl oxidation. Ph.D. thesis, Massachusetts Institute of Technology. (2005). URL: http://hdl.handle.net/1721.1/33716 6. Wächter, A. and Biegler, L.T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming.* 106, 25-57 (2006). DOI: [10.1007/s10107-004-0559-y](https://doi.org/10.1007/s10107-004-0559-y) +7. Pulsipher, J.L., Zhang, W., Hongisto, T.J., and Zavala, V.M. A unifying modeling abstraction for infinite-dimensional optimization. *Computers & Chemical Engineering.* 156, 107567 (2022). DOI: [10.1016/j.compchemeng.2021.107567](https://doi.org/10.1016/j.compchemeng.2021.107567) diff --git a/examples/infinite_model.ipynb b/examples/infinite_model.ipynb new file mode 100644 index 0000000..3153b25 --- /dev/null +++ b/examples/infinite_model.ipynb @@ -0,0 +1,422 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f63f58aa", + "metadata": {}, + "source": [ + "# Kinetic Parameter Estimation\n", + "\n", + "A kinetic parameter estimation problem originally described by [[1](#references)] is used to demonstrate the use of EOptInterface to formulate an `InfiniteOpt.InfiniteModel` [[2](#references)] from an ODE `ModelingToolkit.System` [[3](#references)]. The reacting system of interest is represented by the following initial value problem:\n", + "$$\n", + "\\begin{align*}\n", + " &\\frac{dx_A}{dt} = k_1 x_Z x_Y - c_{\\text{O}_{2}} (k_{2f} + k_{3f}) x_A + \\frac{k_{2f}}{K_2} x_D + \\frac{k_{3f}}{K_3} x_B - k_5 x_A^2 \\\\\n", + " &\\frac{dx_B}{dt} = c_{\\text{O}_{2}} k_{3f} x_A - \\left(\\frac{k_{3f}}{K_3} + k_4\\right) x_B \\\\\n", + " &\\frac{dx_D}{dt} = c_{\\text{O}_{2}} k_{2f} x_A - \\frac{k_{2f}}{K_2} x_D \\\\\n", + " &\\frac{dx_Y}{dt} = -k_1 x_Z x_Y \\\\\n", + " &\\frac{dx_Z}{dt} = -k_1 x_Z x_Y \\\\\n", + " &x_A(0) = x_B(0) = x_D(0) = 0, \\quad x_Y(0) = 0.4, \\quad x_Z(0) = 140,\n", + "\\end{align*}\n", + "$$\n", + "where $x_j$ is the concentration of species $j ∈ \\{A, B, D, Y, Z \\}$ and $T = 273, K_2 = 46\\exp(6500/T-18), K_3 = 2 K_2, k_1 = 53, k_{1s} = k_1 × 10^{-6}, k_5 = 1.2 × 10^{-3},$ and $c_{\\text{O}_{2}} = 2 × 10^{-3}$ are known constants. The remaining uncertain reaction rate constants $\\mathbf{p} = (k_{2f}, k_{3f}, k_4)$ lie within the interval $P = [10,1200] × [10,1200] × [0.001,40]$. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "341eed5e", + "metadata": {}, + "outputs": [], + "source": [ + "using CSV\n", + "using DataFrames\n", + "using EOptInterface\n", + "using InfiniteOpt\n", + "using Ipopt\n", + "using ModelingToolkit\n", + "using ModelingToolkit: t_nounits as t, D_nounits as D" + ] + }, + { + "cell_type": "markdown", + "id": "ac5da559", + "metadata": {}, + "source": [ + "This ODE system is modeled in `ModelingToolkit` with the unknown parameters and initial conditions defined as shown below." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "28356dee", + "metadata": {}, + "outputs": [], + "source": [ + "@mtkmodel KineticParameterEstimation begin\n", + " @parameters begin\n", + " # Known parameters\n", + " T = 273\n", + " K_2 = 46*exp(6500/T-18)\n", + " K_3 = 2*K_2\n", + " k_1 = 53\n", + " k_1s = k_1*1e-6\n", + " k_5 = 1.2e-3\n", + " c_O2 = 2e-3\n", + "\n", + " # Unknown parameters (free design variables)\n", + " k_2f\n", + " k_3f\n", + " k_4\n", + " end\n", + " @variables begin\n", + " # Initial conditions given for differential variables\n", + " x_A(t) = 0.0\n", + " x_B(t) = 0.0\n", + " x_D(t) = 0.0\n", + " x_Y(t) = 0.4\n", + " x_Z(t) = 140.0\n", + " \n", + " I(t)\n", + " end\n", + " @equations begin\n", + " D(x_A) ~ k_1*x_Z*x_Y - c_O2*(k_2f + k_3f)*x_A + k_2f/K_2*x_D + k_3f/K_3*x_B - k_5*x_A^2\n", + " D(x_B) ~ c_O2*k_3f*x_A - (k_3f/K_3 + k_4)*x_B\n", + " D(x_D) ~ c_O2*k_2f*x_A - k_2f/K_2*x_D\n", + " D(x_Y) ~ -k_1s*x_Z*x_Y\n", + " D(x_Z) ~ -k_1*x_Z*x_Y\n", + " I ~ x_A + 2/21*x_B + 2/21*x_D\n", + " end\n", + "end\n", + "\n", + "@mtkcompile system = KineticParameterEstimation();" + ] + }, + { + "cell_type": "markdown", + "id": "a16811dc", + "metadata": {}, + "source": [ + "The objective is to determine the values of uncertain kinetic parameters that best fit this model against experimental intensity versus time data available in `kinetic_intensity_data.csv`, where the relationship between intensity and concentration is known to be $I^{\\text{calc}} = x_A + \\frac{2}{21} x_B + \\frac{2}{21} x_D$ [[4](#references)]. The experimental data set $I^{\\text{exp}}$ contains 201 measurements between $(0, 2)$ spaced linearly by $\\Delta t = 0.01$. Based on this information, we define the integration time span to be $(0, 2)$ to include the initial condition, and use the same time step as between measurements, resulting in $N = 201$ discrete time points. The objective function is then defined as sum of squared error between the model-predicted and measured intensity:\n", + "$$\n", + "\\begin{align}\n", + " f(\\mathbf{x}) = \\sum_{i=2}^{N} \\left(I^{\\text{calc}}(\\mathbf{x}_{i}) - I^{\\text{exp}}_{i} \\right)^2.\n", + "\\end{align}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "413b1871", + "metadata": {}, + "outputs": [], + "source": [ + "data = CSV.read(\"kinetic_intensity_data.csv\", DataFrame)\n", + "intensity(x_A, x_B, x_D) = x_A + 2/21*x_B + 2/21*x_D\n", + "tspan = [data.time[1], data.time[end]]\n", + "N = length(data.time);" + ] + }, + { + "cell_type": "markdown", + "id": "9f776dac", + "metadata": {}, + "source": [ + "We can then define the `InfiniteOpt.InfiniteModel` using an appropriate solver of choice. Given that this is a nonlinear system, we have chosen to use `Ipopt` [[5](#references)]." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "5ca7a46f", + "metadata": {}, + "outputs": [], + "source": [ + "model = InfiniteModel(Ipopt.Optimizer);" + ] + }, + { + "cell_type": "markdown", + "id": "89529b5f", + "metadata": {}, + "source": [ + "We can use `InfiniteOpt.@infinite_parameter` to define our time span and number of supports, while also selecting our desired numerical method to evaluate derivatives." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "c3422460", + "metadata": {}, + "outputs": [], + "source": [ + "@infinite_parameter(model, τ ∈ tspan, num_supports = N, derivative_method = OrthogonalCollocation(2));" + ] + }, + { + "cell_type": "markdown", + "id": "354befae", + "metadata": {}, + "source": [ + "We can use `EOptInterface.decision_vars` on the `ModelingToolkit.System` to retrieve the decision variables for the optimization problem." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "8f38e3b9", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "8-element Vector{SymbolicUtils.BasicSymbolicImpl.var\"typeof(BasicSymbolicImpl)\"{SymReal}}:\n", + " x_Z(t)\n", + " x_Y(t)\n", + " x_D(t)\n", + " x_B(t)\n", + " x_A(t)\n", + " k_2f\n", + " k_3f\n", + " k_4" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "decision_vars(system)" + ] + }, + { + "cell_type": "markdown", + "id": "69d397b1", + "metadata": {}, + "source": [ + "Now that we know which variables appear in the system, we can add the decision variables to the `InfiniteOpt.InfiniteModel` with appropriate bounds." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "7a7f6375", + "metadata": {}, + "outputs": [], + "source": [ + "V = length(unknowns(system))\n", + "@variable(model, -75.0 <= z[1:V] <= 150.0, Infinite(τ))\n", + "pL = [10.0, 10.0, 0.001]\n", + "pU = [1200.0, 1200.0, 40.0]\n", + "@variable(model, pL[i] <= p[i=1:3] <= pU[i]);" + ] + }, + { + "cell_type": "markdown", + "id": "fdeb4018", + "metadata": {}, + "source": [ + "We then use `EOptInterface.register_odesystem` to register the ODE `ModelingToolkit.System` using explicit Euler and register the system as differential constraints." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "9b15d798", + "metadata": {}, + "outputs": [], + "source": [ + "register_odesystem(model, system)" + ] + }, + { + "cell_type": "markdown", + "id": "cef6f20c", + "metadata": {}, + "source": [ + "We can add the objective function defined in Equation 1 to the `InfiniteOpt.InfiniteModel`." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "d9762428", + "metadata": {}, + "outputs": [], + "source": [ + "@objective(model, Min, sum((intensity(z[5](i), z[4](i), z[3](i)) - data.intensity[findfirst(==(i), data.time)])^2 for i in supports(τ)));" + ] + }, + { + "cell_type": "markdown", + "id": "5d481ed3", + "metadata": {}, + "source": [ + "Finally, we optimize the `InfiniteOpt.InfiniteModel` and retrieve the results." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "a04a5a8e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\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 https://github.com/coin-or/Ipopt\n", + "******************************************************************************\n", + "\n", + "This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.\n", + "\n", + "Number of nonzeros in equality constraint Jacobian...: 9035\n", + "Number of nonzeros in inequality constraint Jacobian.: 0\n", + "Number of nonzeros in Lagrangian Hessian.............: 3819\n", + "\n", + "Total number of variables............................: 2013\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 1008\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 2010\n", + "Total number of inequality constraints...............: 0\n", + " inequality constraints with only lower bounds: 0\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 4.7973855e+05 1.40e+02 1.00e+02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 4.7973855e+05 7.42e+02 1.18e+05 -1.0 1.40e+02 - 3.46e-01 5.00e-01f 2\n", + " 2 4.7653605e+05 7.41e+02 1.18e+05 -1.0 2.23e+03 - 2.65e-01 1.21e-03f 1\n", + " 3 3.3502461e+05 5.92e+02 8.94e+04 -1.0 2.23e+03 - 2.76e-01 2.80e-01f 1\n", + " 4 2.0021793e+06 3.78e+02 4.13e+04 -1.0 1.96e+03 - 6.41e-02 9.59e-01h 1\n", + " 5 1.4085380e+06 5.17e+00 4.31e+03 -1.0 4.64e+02 - 7.76e-03 1.00e+00f 1\n", + " 6 1.3484374e+06 1.43e+00 2.64e+03 -1.0 4.83e+00 - 1.87e-02 1.00e+00f 1\n", + " 7 1.3095850e+06 9.59e-01 2.58e+03 -1.0 3.61e+02 - 1.33e-02 2.60e-02f 1\n", + " 8 1.1557650e+06 2.51e+00 2.27e+03 -1.0 2.06e+02 - 7.04e-02 1.21e-01f 1\n", + " 9 3.2268278e+05 6.40e+01 4.15e+02 -1.0 2.02e+02 - 5.05e-01 1.00e+00f 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 7.8612401e+04 1.45e+01 9.31e+01 -1.0 2.66e+02 - 5.66e-01 1.00e+00f 1\n", + " 11 7.0581018e+04 1.74e+00 8.74e+01 -1.0 8.23e+01 0.0 2.58e-01 1.00e+00f 1\n", + " 12 4.8543183e+04 1.45e+00 5.38e+01 -1.0 1.61e+02 - 9.39e-01 4.48e-01f 1\n", + " 13 2.5469240e+04 3.41e+00 1.39e+01 -1.0 1.62e+02 - 7.62e-01 1.00e+00f 1\n", + " 14 2.1293158e+04 1.50e+00 4.61e+00 -1.0 1.07e+02 - 6.61e-03 1.00e+00f 1\n", + " 15 1.6956787e+04 3.38e+01 2.91e+00 -1.0 1.03e+03 - 2.75e-02 1.00e+00f 1\n", + " 16 1.6838348e+04 2.36e+00 1.05e+00 -1.0 2.25e+02 - 6.16e-01 1.00e+00f 1\n", + " 17 1.6831934e+04 9.40e-02 4.83e-01 -1.0 4.16e+01 - 5.34e-01 1.00e+00h 1\n", + " 18 1.6831563e+04 1.63e-02 1.77e-01 -1.0 5.47e-01 -0.5 1.00e+00 1.00e+00h 1\n", + " 19 1.6826536e+04 9.12e-01 4.91e-01 -2.5 4.12e+00 - 2.33e-01 1.00e+00f 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 20 1.6825176e+04 6.33e-02 7.73e-02 -2.5 8.06e-01 -1.0 1.00e+00 1.00e+00h 1\n", + " 21 1.6779023e+04 7.95e+01 2.34e+00 -2.5 4.95e+01 - 1.24e-01 6.81e-01f 1\n", + " 22 1.6778660e+04 7.87e+01 2.32e+00 -2.5 3.79e+02 - 8.23e-02 9.19e-03H 1\n", + " 23 1.6809951e+04 6.85e+01 1.44e+01 -2.5 5.79e+01 -1.4 8.14e-05 1.25e-01h 4\n", + " 24 1.6818131e+04 6.64e+01 1.49e+01 -2.5 5.30e+01 -1.9 1.00e+00 3.12e-02h 6\n", + " 25 1.6824219e+04 6.43e+01 1.48e+01 -2.5 2.31e+01 -1.5 1.33e-01 3.12e-02h 6\n", + " 26 1.6825654e+04 6.38e+01 1.48e+01 -2.5 3.07e+01 -1.1 1.00e+00 7.81e-03h 8\n", + " 27 1.6826248e+04 6.36e+01 1.47e+01 -2.5 2.74e+01 -1.5 1.36e-01 3.91e-03h 9\n", + " 28 1.6757369e+04 1.26e+02 1.20e+03 -2.5 2.72e+01 -1.1 1.00e+00 1.00e+00H 1\n", + " 29 1.6804218e+04 1.21e+02 1.16e+03 -2.5 3.07e+03 - 1.45e-01 3.02e-02h 4\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 30 1.6823590e+04 1.12e+02 1.07e+03 -2.5 6.25e+01 -1.6 1.00e+00 7.84e-02h 4\n", + " 31 1.6825121e+04 1.11e+02 1.07e+03 -2.5 1.91e+02 -2.1 1.00e+00 3.49e-03h 8\n", + " 32 1.6825468e+04 1.11e+02 1.07e+03 -2.5 6.58e+02 -2.5 1.00e+00 2.96e-04h 11\n", + " 33 1.6825622e+04 1.11e+02 1.07e+03 -2.5 3.16e+03 -3.0 5.62e-01 3.35e-05h 13\n", + " 34 1.6825715e+04 1.11e+02 1.07e+03 -2.5 7.72e+02 -2.6 1.00e+00 6.98e-05h 13\n", + " 35 1.6934202e+04 1.09e+02 1.05e+03 -2.5 3.85e+03 -3.1 1.64e-01 1.49e-02h 4\n", + " 36 1.7117832e+04 9.91e+01 1.01e+03 -2.5 8.85e+02 -2.6 1.00e+00 6.57e-02h 3\n", + " 37 1.8844669e+04 4.87e+01 3.21e+03 -2.5 1.77e+03 - 1.00e+00 1.00e+00H 1\n", + " 38 1.8000776e+04 6.96e+00 3.19e+03 -2.5 1.54e+02 -3.1 6.49e-01 1.00e+00f 1\n", + " 39 1.7541612e+04 1.85e+00 2.82e+02 -2.5 2.60e+01 -1.8 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", + " 40 1.7095381e+04 5.12e+00 1.41e+02 -2.5 1.14e+03 - 3.82e-03 5.00e-01f 2\n", + " 41 1.7017993e+04 1.83e+00 1.97e+01 -2.5 9.39e+01 - 1.00e+00 1.00e+00f 1\n", + " 42 1.7015195e+04 5.78e-02 3.74e-01 -2.5 9.86e-01 -2.3 1.00e+00 1.00e+00h 1\n", + " 43 1.7012013e+04 9.52e-02 1.82e-01 -2.5 6.85e+00 -2.7 1.00e+00 1.00e+00h 1\n", + " 44 1.7000452e+04 9.70e-01 6.59e-01 -2.5 3.24e+01 -3.2 1.00e+00 1.00e+00f 1\n", + " 45 1.6989569e+04 4.13e-01 1.05e-01 -2.5 1.41e+01 -2.8 1.00e+00 1.00e+00h 1\n", + " 46 1.6932155e+04 2.02e+01 8.48e+00 -2.5 9.97e+01 -3.3 1.00e+00 1.00e+00f 1\n", + " 47 1.6816745e+04 1.71e+01 9.57e+00 -2.5 1.89e+01 -2.8 1.00e+00 1.00e+00f 1\n", + " 48 1.6826444e+04 6.15e+00 3.95e+00 -2.5 7.61e+00 - 1.00e+00 1.00e+00h 1\n", + " 49 1.6774478e+04 1.93e+01 2.48e-01 -2.5 1.78e+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", + " 50 1.6797859e+04 3.81e+00 3.67e-01 -2.5 6.83e+00 - 1.00e+00 1.00e+00h 1\n", + " 51 1.6795905e+04 2.88e-01 5.15e-03 -2.5 2.44e+00 - 1.00e+00 1.00e+00h 1\n", + " 52 1.6796030e+04 3.96e-03 8.31e-05 -2.5 2.18e-01 - 1.00e+00 1.00e+00h 1\n", + " 53 1.6796032e+04 3.33e-07 1.46e-08 -3.8 4.11e-02 - 1.00e+00 1.00e+00h 1\n", + " 54 1.6796032e+04 1.61e-09 8.79e-11 -5.7 2.29e-03 - 1.00e+00 1.00e+00h 1\n", + " 55 1.6796032e+04 3.53e-13 1.93e-12 -8.6 2.83e-05 - 1.00e+00 1.00e+00h 1\n", + "\n", + "Number of Iterations....: 55\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: 6.6954341638766082e+03 1.6796032234817583e+04\n", + "Dual infeasibility......: 1.9321875863753234e-12 4.8470471354294094e-12\n", + "Constraint violation....: 3.5260683262094972e-13 3.5260683262094972e-13\n", + "Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Complementarity.........: 2.5060333824815601e-09 6.2865852226255921e-09\n", + "Overall NLP error.......: 2.5060333824815601e-09 6.2865852226255921e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 171\n", + "Number of objective gradient evaluations = 56\n", + "Number of equality constraint evaluations = 171\n", + "Number of inequality constraint evaluations = 0\n", + "Number of equality constraint Jacobian evaluations = 56\n", + "Number of inequality constraint Jacobian evaluations = 0\n", + "Number of Lagrangian Hessian evaluations = 55\n", + "Total seconds in IPOPT = 0.379\n", + "\n", + "EXIT: Optimal Solution Found.\n", + "Termination Status: LOCALLY_SOLVED\n", + "Primal Status: FEASIBLE_POINT\n", + "Solve Time: 0.388\n", + "f^* = 16796.03223\n", + "p* = [798.186, 423.947, 12.977]\n" + ] + } + ], + "source": [ + "optimize!(model)\n", + "println(\"Termination Status: $(termination_status(model))\")\n", + "println(\"Primal Status: $(primal_status(model))\")\n", + "println(\"Solve Time: $(round.(solve_time(model), digits=5))\")\n", + "println(\"f^* = $(round(objective_value(model), digits=5))\")\n", + "println(\"p* = $(round.(value.(p), digits=3))\")" + ] + }, + { + "cell_type": "markdown", + "id": "4ef1459e", + "metadata": {}, + "source": [ + "### References\n", + "\n", + "1. Taylor, J.W. Direct measurement and analysis of cyclohexadienyl oxidation. Ph.D. thesis, Massachusetts Institute of Technology. (2005). URL: http://hdl.handle.net/1721.1/33716\n", + "2. Pulsipher, J.L., Zhang, W., Hongisto, T.J., and Zavala, V.M. A unifying modeling abstraction for infinite-dimensional optimization. *Computers & Chemical Engineering.* 156, 107567 (2022). DOI: [10.1016/j.compchemeng.2021.107567](https://doi.org/10.1016/j.compchemeng.2021.107567)\n", + "3. Ma, Y., Gowda, S., Anantharaman, R., Laughman, C., Shah, V., and Rackauckas, C. ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling. (2022). DOI: [10.48550/arXiv.2103.05244](https://doi.org/10.48550/arXiv.2103.05244)\n", + "3. Singer, A.B. Global dynamic optimization. Ph.D. thesis, Massachusetts Institute of Technology. (2004). URL: http://hdl.handle.net/1721.1/28662\n", + "4. Wächter, A. and Biegler, L.T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming.* 106, 25-57 (2006). DOI: [10.1007/s10107-004-0559-y](https://doi.org/10.1007/s10107-004-0559-y)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.11", + "language": "julia", + "name": "julia-1.10" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/infinite_model.jl b/examples/infinite_model.jl new file mode 100644 index 0000000..ba4f8d3 --- /dev/null +++ b/examples/infinite_model.jl @@ -0,0 +1,89 @@ +using CSV +using DataFrames +using EOptInterface +using InfiniteOpt +using Ipopt +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D + +# Create ModelingToolkit ODE system +@mtkmodel KineticParameterEstimation begin + @parameters begin + # Known parameters + T = 273.0 + K_2 = 46.0*exp(6500.0/T - 18.0) + K_3 = 2.0*K_2 + k_1 = 53.0 + k_1s = k_1*1e-6 + k_5 = 1.2e-3 + c_O2 = 2e-3 + + # Unknown parameters (free design variables) + k_2f + k_3f + k_4 + end + @variables begin + # Initial conditions given for differential variables + x_A(t) = 0.0 + x_B(t) = 0.0 + x_D(t) = 0.0 + x_Y(t) = 0.4 + x_Z(t) = 140.0 + + I(t) + end + @equations begin + D(x_A) ~ k_1*x_Z*x_Y - c_O2*(k_2f + k_3f)*x_A + k_2f/K_2*x_D + k_3f/K_3*x_B - k_5*x_A^2 + D(x_B) ~ c_O2*k_3f*x_A - (k_3f/K_3 + k_4)*x_B + D(x_D) ~ c_O2*k_2f*x_A - k_2f/K_2*x_D + D(x_Y) ~ -k_1s*x_Z*x_Y + D(x_Z) ~ -k_1*x_Z*x_Y + I ~ x_A + 2/21*x_B + 2/21*x_D + end +end + +# Compile system +@mtkcompile system = KineticParameterEstimation() + +# Include experimental intensity data +data = CSV.read("examples/kinetic_intensity_data.csv", DataFrame) +# Define intensity function +intensity(x_A, x_B, x_D) = x_A + 2/21*x_B + 2/21*x_D +# Retrieve integration time data +tspan = [data.time[1], data.time[end]] +tstep = data.time[2] - data.time[1] +N = length(data.time) + +# Create InfiniteOpt model +model = InfiniteModel(Ipopt.Optimizer) + +# Retrieve decision variables from ModelingToolkit system +# Returns [x_Z(t), x_Y(t), x_D(t), x_B(t), x_A(t), k_2f, k_3f, k_4] +decision_vars(system) + +# Create state decision variables +# z = [x_Z(t), x_Y(t), x_D(t), x_B(t), x_A(t)] +V = length(unknowns(system)) +@infinite_parameter(model, τ ∈ tspan, num_supports = N, derivative_method = OrthogonalCollocation(2)) +@variable(model, -75.0 <= z[1:V] <= 150.0, Infinite(τ)) + +# Create free design decision variables +# p = [k_2f, k_3f, k_4] +pL = [10.0, 10.0, 0.001] +pU = [1200.0, 1200.0, 40.0] +@variable(model, pL[i] <= p[i=1:3] <= pU[i]) + +# Register ModelingToolkit ODE system as constraints +register_odesystem(model, system) + +# Define objective function +@objective(model, Min, sum((intensity(z[5](i), z[4](i), z[3](i)) - data.intensity[findfirst(==(i), data.time)])^2 for i in supports(τ))) + +# Optimize model and retrieve results +optimize!(model) +println("Termination Status: $(termination_status(model))") +println("Primal Status: $(primal_status(model))") +println("Solve Time: $(round.(solve_time(model), digits=5))") +println("f^* = $(round(objective_value(model), digits=5))") +println("p* = $(round.(value.(p), digits=3))") diff --git a/ext/EOIInfiniteOptExt.jl b/ext/EOIInfiniteOptExt.jl new file mode 100644 index 0000000..5314f85 --- /dev/null +++ b/ext/EOIInfiniteOptExt.jl @@ -0,0 +1,65 @@ +# Copyright (c) 2025 Joseph Choi, Dimitri Alston, Pengfei Xu, Matthew Stuber, +# and the University of Connecticut (UConn) +# This code is licensed under the MIT license (see LICENSE.md for full details). +################################################################################ +# EOptInterface +# An abstraction layer for optimizing equation-oriented/acausal models +# https://github.com/PSORLab/EOptInterface.jl +################################################################################ +# ext/EOIInfiniteOptExt.jl +# An extension for using InfiniteOpt. +################################################################################ + +module EOIInfiniteOptExt + + # Imports + import DocStringExtensions + import EOptInterface + import InfiniteOpt + import JuMP + import ModelingToolkit + import SymbolicUtils + import Symbolics + + # Functions + """ + $(DocStringExtensions.TYPEDSIGNATURES) + + Automatically registers ODEs from a `ModelingToolkit.System` as + InfiniteOpt constraints. + """ + function EOptInterface.register_odesystem(model::InfiniteOpt.InfiniteModel, sys::ModelingToolkit.System) + t = InfiniteOpt.all_parameters(model, InfiniteOpt.InfiniteParameter)[1] + # Number of ODE variables + V = length(ModelingToolkit.unknowns(sys)) + param_dict = copy(ModelingToolkit.initial_conditions(sys).dict) + for var in ModelingToolkit.unknowns(sys) + pop!(param_dict, var) + end + dx = [] + for j in 1:V + dxj_expr = ModelingToolkit.full_equations(sys)[j].rhs + dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(sys)) + dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(sys)) + # Fully substitute parameters with default values + while ~isempty(intersect(Symbolics.get_variables(dxj_expr), keys(param_dict))) + dxj_expr = SymbolicUtils.substitute(dxj_expr, param_dict) + end + dxj = Symbolics.build_function( + dxj_expr, + EOptInterface.decision_vars(sys)..., + expression = Val{false} + ) + push!(dx, dxj) + end + ps = JuMP.all_variables(model)[end-length(setdiff(EOptInterface.decision_vars(sys), ModelingToolkit.unknowns(sys)))+1:end] + xs = reshape(setdiff(JuMP.all_variables(model), ps), V) + # Extract initial conditions from the ModelingToolkit system + ic_dict = [ModelingToolkit.initial_conditions(sys)[i].val for i in intersect(ModelingToolkit.unknowns(sys), keys(ModelingToolkit.initial_conditions(sys)))] + # Formulate InfiniteOpt constraints + InfiniteOpt.@constraint(model, [i in 1:V], xs[i](0) == ic_dict[i]) + InfiniteOpt.@constraint(model, [i in 1:V], InfiniteOpt.∂(xs[i], t) == dx[i](xs..., ps...)) + return + end + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 8e19776..1ff5395 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using CSV using DataFrames using EOptInterface +using InfiniteOpt using Ipopt using JuMP using ModelingToolkit @@ -232,3 +233,62 @@ end @test isapprox(JuMP.objective_value(model), 16796.032234817612, atol=1e-3) end + +@testset "InfiniteOpt ODE Model" begin + + @mtkmodel KineticParameterEstimation begin + @parameters begin + T = 273.0 + K_2 = 46.0*exp(6500.0/T - 18.0) + K_3 = 2.0*K_2 + k_1 = 53.0 + k_1s = k_1*1e-6 + k_5 = 1.2e-3 + c_O2 = 2e-3 + k_2f + k_3f + k_4 + end + @variables begin + x_A(t) = 0.0 + x_B(t) = 0.0 + x_D(t) = 0.0 + x_Y(t) = 0.4 + x_Z(t) = 140.0 + I(t) + end + @equations begin + D(x_A) ~ k_1*x_Z*x_Y - c_O2*(k_2f + k_3f)*x_A + k_2f/K_2*x_D + k_3f/K_3*x_B - k_5*x_A^2 + D(x_B) ~ c_O2*k_3f*x_A - (k_3f/K_3 + k_4)*x_B + D(x_D) ~ c_O2*k_2f*x_A - k_2f/K_2*x_D + D(x_Y) ~ -k_1s*x_Z*x_Y + D(x_Z) ~ -k_1*x_Z*x_Y + I ~ x_A + 2/21*x_B + 2/21*x_D + end + end + + @mtkcompile system = KineticParameterEstimation() + + data = CSV.read("kinetic_intensity_data.csv", DataFrame) + intensity(x_A, x_B, x_D) = x_A + 2/21*x_B + 2/21*x_D + tspan = [data.time[1], data.time[end]] + N = length(data.time) + + model = InfiniteOpt.InfiniteModel(Ipopt.Optimizer) + V = length(unknowns(system)) + InfiniteOpt.@infinite_parameter(model, τ ∈ tspan, num_supports = N, derivative_method = OrthogonalCollocation(2)) + InfiniteOpt.@variable(model, -75.0 <= z[1:V] <= 150.0, Infinite(τ)) + pL = [10.0, 10.0, 0.001] + pU = [1200.0, 1200.0, 40.0] + JuMP.@variable(model, pL[i] <= p[i=1:3] <= pU[i]) + + EOptInterface.register_odesystem(model, system) + + JuMP.@objective(model, Min, sum((intensity(z[5](i), z[4](i), z[3](i)) - data.intensity[findfirst(==(i), data.time)])^2 for i in supports(τ))) + JuMP.optimize!(model) + + @test JuMP.termination_status(model) == JuMP.LOCALLY_SOLVED + @test JuMP.primal_status(model) == JuMP.FEASIBLE_POINT + @test isapprox(JuMP.objective_value(model), 16796.032234817612, atol=1e-3) + +end From 4bdcb8927af6dfb00e675eb5e04edc21a52d6c24 Mon Sep 17 00:00:00 2001 From: Dimitri Alston <123396563+DimitriAlston@users.noreply.github.com> Date: Thu, 25 Jun 2026 11:56:48 -0400 Subject: [PATCH 2/3] Pass infinite parameter as input to `register_odesystem` --- examples/infinite_model.ipynb | 4 ++-- examples/infinite_model.jl | 2 +- ext/EOIInfiniteOptExt.jl | 3 +-- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/infinite_model.ipynb b/examples/infinite_model.ipynb index 3153b25..9ce1c42 100644 --- a/examples/infinite_model.ipynb +++ b/examples/infinite_model.ipynb @@ -221,12 +221,12 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "id": "9b15d798", "metadata": {}, "outputs": [], "source": [ - "register_odesystem(model, system)" + "register_odesystem(model, system, τ)" ] }, { diff --git a/examples/infinite_model.jl b/examples/infinite_model.jl index ba4f8d3..fb9f04b 100644 --- a/examples/infinite_model.jl +++ b/examples/infinite_model.jl @@ -75,7 +75,7 @@ pU = [1200.0, 1200.0, 40.0] @variable(model, pL[i] <= p[i=1:3] <= pU[i]) # Register ModelingToolkit ODE system as constraints -register_odesystem(model, system) +register_odesystem(model, system, τ) # Define objective function @objective(model, Min, sum((intensity(z[5](i), z[4](i), z[3](i)) - data.intensity[findfirst(==(i), data.time)])^2 for i in supports(τ))) diff --git a/ext/EOIInfiniteOptExt.jl b/ext/EOIInfiniteOptExt.jl index 5314f85..725be47 100644 --- a/ext/EOIInfiniteOptExt.jl +++ b/ext/EOIInfiniteOptExt.jl @@ -28,8 +28,7 @@ module EOIInfiniteOptExt Automatically registers ODEs from a `ModelingToolkit.System` as InfiniteOpt constraints. """ - function EOptInterface.register_odesystem(model::InfiniteOpt.InfiniteModel, sys::ModelingToolkit.System) - t = InfiniteOpt.all_parameters(model, InfiniteOpt.InfiniteParameter)[1] + function EOptInterface.register_odesystem(model::InfiniteOpt.InfiniteModel, sys::ModelingToolkit.System, t::InfiniteOpt.GeneralVariableRef) # Number of ODE variables V = length(ModelingToolkit.unknowns(sys)) param_dict = copy(ModelingToolkit.initial_conditions(sys).dict) From 085d8e97da38b7848ce2cf3ea03dbd91a4cdc21e Mon Sep 17 00:00:00 2001 From: Dimitri Alston <123396563+DimitriAlston@users.noreply.github.com> Date: Thu, 25 Jun 2026 12:11:18 -0400 Subject: [PATCH 3/3] Update runtests.jl --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1ff5395..0219f8d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -282,7 +282,7 @@ end pU = [1200.0, 1200.0, 40.0] JuMP.@variable(model, pL[i] <= p[i=1:3] <= pU[i]) - EOptInterface.register_odesystem(model, system) + EOptInterface.register_odesystem(model, system, τ) JuMP.@objective(model, Min, sum((intensity(z[5](i), z[4](i), z[3](i)) - data.intensity[findfirst(==(i), data.time)])^2 for i in supports(τ))) JuMP.optimize!(model)