Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ authors = ["Joseph Choi <jsphchoi@mit.edu>", "Dimitri Alston <dimitri.alston@uco
version = "0.2.0"

[deps]
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
Expand All @@ -12,6 +13,7 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

[compat]
DocStringExtensions = "0.9"
JuMP = "1"
ModelingToolkit = "11"
Reexport = "1"
Expand Down
27 changes: 2 additions & 25 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,40 +10,17 @@
|:-----------------------------------------------------------------:|
| [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://PSORLab.github.io/EOptInterface.jl/stable) [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://PSORLab.github.io/EOptInterface.jl/dev) |

## Feature Summary

```julia
decision_vars(::ModelingToolkit.System)
```
Returns the decision variables for an optimization problem formulated from a `ModelingToolkit` system.

```julia
register_nlsystem(::JuMP.Model, ::ModelingToolkit.System, obj::Symbolics.Num, ineqs::Vector{Symbolics.Num})
```
Automatically formulates algebraic `JuMP` constraints and objective function from
an algebraic `ModelingToolkit` system and user-provided constraints and objective symbolic expressions.

```julia
full_solution(::JuMP.Model, ::ModelingToolkit.System)
```
Returns a dictionary of optimal solution values for the observed variables of an algebraic `ModelingToolkit` system if the `JuMP` model is solved.

```julia
register_odesystem(::JuMP.Model, ::ModelingToolkit.System, tspan::Tuple{Real,Real}, tstep::Real, solver::String)
```
Automatically applies forward transcription and registers the discretized ODE `ModelingToolkit` system as algebraic `JuMP` constraints. Available integration methods: `"EE"` (explicit Euler), `"IE"` (implicit Euler).

## Examples

The code for these examples can be found in the [`examples/`](https://github.com/PSORLab/EOptInterface.jl/tree/main/examples) subdirectory.
The code and notebooks for these examples can be found in the [`examples/`](https://github.com/PSORLab/EOptInterface.jl/tree/main/examples) subdirectory.

### [Algebraic System](https://github.com/PSORLab/EOptInterface.jl/blob/main/examples/algebraic_model.jl)

An optimal reactor-separator-recycle process design problem originally presented by [[3](#references)] is used to demonstrate the use of `register_nlsystem` to formulate and solve a reduced-space model using the deterministic global optimizer [`EAGO.jl`](https://github.com/PSORLab/EAGO.jl) [[4](#references)].

### [ODE System](https://github.com/PSORLab/EOptInterface.jl/blob/main/examples/ode_model.jl)

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` [[6](#references)].
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)].

## 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)
Expand Down
3 changes: 2 additions & 1 deletion src/EOptInterface.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,13 @@
module EOptInterface

# Imports
import DocStringExtensions
import JuMP
import ModelingToolkit
import ModelingToolkit: t_nounits as t, D_nounits as D
import Reexport
import Symbolics
import SymbolicUtils
import Symbolics

# Reexports
Reexport.@reexport using SciCompDSL
Expand Down
56 changes: 24 additions & 32 deletions src/userfuncs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,10 @@
################################################################################

"""
decision_vars(sys)
$(DocStringExtensions.TYPEDSIGNATURES)

Returns the decision variables for an optimization problem formulated from a ModelingToolkit system.
Returns the decision variables for an optimization problem from a
`ModelingToolkit.System`.
"""
function decision_vars(sys::ModelingToolkit.System)
return [
Expand All @@ -25,16 +26,11 @@ function decision_vars(sys::ModelingToolkit.System)
end

"""
register_nlsystem(model, sys, obj, ineqs)
$(DocStringExtensions.TYPEDSIGNATURES)

Automatically formulates algebraic JuMP constraints and objective function from
an algebraic ModelingToolkit system and user-provided constraints and objective symbolic expressions.

# Arguments
- `model::JuMP.Model`: the JuMP model
- `sys::ModelingToolkit.System`: the ModelingToolkit system
- `obj::Symbolics.Num`: a symbolic expression of the objective function using the ModelingToolkit system variables
- `ineqs::Vector{Symbolics.Num}`: a vector of symbolic expressions of inequality constraints using the ModelingToolkit system variables
Automatically formulates and adds user-provided `Symbolics.Num` objective
function and `Vector{Symbolics.Num}` constraints from a `ModelingToolkit.System`
to a `JuMP.Model`.
"""
function register_nlsystem(model::JuMP.Model, sys::ModelingToolkit.System, obj::Symbolics.Num, ineqs::Vector{Symbolics.Num})
h = EOptInterface.mtk_generate_model_equations(sys)
Expand All @@ -51,49 +47,44 @@ function register_nlsystem(model::JuMP.Model, sys::ModelingToolkit.System, obj::
end

"""
register_odesystem(model, sys, tspan, tstep, integrator)

Automatically applies forward transcription and registers the discretized ODE ModelingToolkit system as algebraic JuMP constraints.
$(DocStringExtensions.TYPEDSIGNATURES)

# Arguments
- `model::JuMP.Model`: the JuMP model
- `sys::ModelingToolkit.System`: the ModelingToolkit system
- `tspan::Tuple{Real,Real}`: the time span over which the dynamic model is simulated
- `tstep::Real`: the time step used in the integration scheme
- `integrator::String`: integration scheme used in discretization, `"EE"` for explicit Euler or `"IE"` for implicit Euler
Automatically applies specified direct transcription method and registers the
discretized ODE `ModelingToolkit.System` as algebraic JuMP constraints.
Current supports Explicit Euler "EE" and Implicit Euler "IE".
"""
function register_odesystem(model::JuMP.Model, odesys::ModelingToolkit.System, tspan::Tuple{Real,Real}, tstep::Real, integrator::String)
function register_odesystem(model::JuMP.Model, sys::ModelingToolkit.System, tspan::Tuple{Real,Real}, tstep::Real, integrator::String)
if integrator != "EE" && integrator != "IE"
error("Available integrators: EE, IE")
end
# Number of discrete time nodes
N = Int(floor((tspan[2] - tspan[1])/tstep)) + 1
# Number of ODE variables
V = length(ModelingToolkit.unknowns(odesys))
param_dict = copy(ModelingToolkit.initial_conditions(odesys).dict)
for var in ModelingToolkit.unknowns(odesys)
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(odesys)[j].rhs
dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(odesys))
dxj_expr = SymbolicUtils.substitute(dxj_expr, ModelingToolkit.bindings(odesys))
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(odesys)...,
EOptInterface.decision_vars(sys)...,
expression = Val{false}
)
push!(dx, dxj)
end
ps = JuMP.all_variables(model)[end-length(setdiff(EOptInterface.decision_vars(odesys), ModelingToolkit.unknowns(odesys)))+1: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, N)
# Extract initial conditions from the ModelingToolkit system and fix them in the JuMP model for x[1:V,1]
JuMP.fix.(xs[:,1], [ModelingToolkit.initial_conditions(odesys)[ModelingToolkit.unknowns(odesys)[i]].val for i in eachindex(ModelingToolkit.unknowns(odesys))], force=true)
JuMP.fix.(xs[:,1], [ModelingToolkit.initial_conditions(sys)[ModelingToolkit.unknowns(sys)[i]].val for i in eachindex(ModelingToolkit.unknowns(sys))], force=true)
# Formulate JuMP constraints based on chosen ODE discretization method
if integrator == "EE"
JuMP.@constraint(model, [j in 1:V, i in 1:(N-1)], xs[j,i+1] == xs[j,i] + tstep*dx[j](xs[:,i]..., ps...))
Expand All @@ -104,9 +95,10 @@ function register_odesystem(model::JuMP.Model, odesys::ModelingToolkit.System, t
end

"""
full_solution(model, sys)
$(DocStringExtensions.TYPEDSIGNATURES)

Returns a dictionary of optimal solution values for the observed variables of an algebraic ModelingToolkit system if the JuMP model is solved.
Returns a dictionary of optimal solution values for the observed variables of a
`ModelingToolkit.System` if the `JuMP.Model` is solved.
"""
function full_solution(model::JuMP.Model, sys::ModelingToolkit.System)
vars = EOptInterface.decision_vars(sys)
Expand Down
Loading