Merge branch 'master' into myb_fb/simplify_clock

Author: YingboMa

.github/workflows/Downstream.yml

@@ -22,7 +22,7 @@ jobs:
           - {user: SciML, repo: CellMLToolkit.jl, group: All}
           - {user: SciML, repo: SBMLToolkit.jl, group: All}
           - {user: SciML, repo: NeuralPDE.jl, group: NNPDE}
-          - {user: SciML, repo: DataDrivenDiffEq.jl, group: Standard}
+          #- {user: SciML, repo: DataDrivenDiffEq.jl, group: Standard}
           - {user: SciML, repo: StructuralIdentifiability.jl, group: All}
           - {user: SciML, repo: ModelingToolkitStandardLibrary.jl}
           - {user: SciML, repo: ModelOrderReduction.jl, group: All}

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "8.30.0"
+version = "8.32.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -72,7 +72,7 @@ NonlinearSolve = "0.3.8"
 RecursiveArrayTools = "2.3"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.4.3, 0.5"
-SciMLBase = "1.60.0"
+SciMLBase = "1.70.0"
 Setfield = "0.7, 0.8, 1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"

README.md

@@ -3,7 +3,7 @@
 
 [![Join the chat at https://julialang.zulipchat.com #sciml-bridged](https://img.shields.io/static/v1?label=Zulip&message=chat&color=9558b2&labelColor=389826)](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged)
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](http://mtk.sciml.ai/stable/)
-[![Global Docs](https://img.shields.io/badge/docs-SciML-blue.svg)](https://docs.sciml.ai/dev/modules/ModelingToolkit/)
+[![Global Docs](https://img.shields.io/badge/docs-SciML-blue.svg)](https://docs.sciml.ai/ModelingToolkit/stable/)
 
 [![codecov](https://codecov.io/gh/SciML/ModelingToolkit.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/SciML/ModelingToolkit.jl)
 [![Coverage Status](https://coveralls.io/repos/github/SciML/ModelingToolkit.jl/badge.svg?branch=master)](https://coveralls.io/github/SciML/ModelingToolkit.jl?branch=master)

docs/Project.toml

@@ -3,11 +3,12 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

docs/src/systems/OptimizationSystem.md

@@ -8,9 +8,10 @@ OptimizationSystem
 
 ## Composition and Accessor Functions
 
-- `get_eqs(sys)` or `equations(sys)`: The equation to be minimized.
+- `get_op(sys)` or `objective(sys)`: The objective to be minimized.
 - `get_states(sys)` or `states(sys)`: The set of states for the optimization.
 - `get_ps(sys)` or `parameters(sys)`: The parameters for the optimization.
+- `get_constraints(sys)` or `constraints(sys)`: The constraints for the optimization.
 
 ## Transformations
 

docs/src/tutorials/optimization.md

@@ -1,26 +1,91 @@
 # Modeling Optimization Problems
 
-```julia
-using ModelingToolkit, Optimization, OptimizationOptimJL
+## Rosenbrock Function in 2D
+Let's solve the classical _Rosenbrock function_ in two dimensions.
 
-@variables x y
-@parameters a b
+First, we need to make some imports.
+```@example rosenbrock_2d
+using ModelingToolkit, Optimization, OptimizationOptimJL
+```
+Now we can define our optimization problem.
+```@example rosenbrock_2d
+@variables begin
+    x, [bounds = (-2.0, 2.0)]
+    y, [bounds = (-1.0, 3.0)]
+end
+@parameters a = 1 b = 1
 loss = (a - x)^2 + b * (y - x^2)^2
-@named sys = OptimizationSystem(loss,[x,y],[a,b])
+@named sys = OptimizationSystem(loss, [x, y], [a, b])
+```
+
+A visualization of the objective function is depicted below.
+```@eval
+using Plots
+x = -2:0.01:2
+y = -1:0.01:3
+contour(x, y, (x,y) -> (1 - x)^2 + 100 * (y - x^2)^2, fill=true, color=:viridis, ratio=:equal, xlims=(-2, 2))
+savefig("obj_fun.png"); nothing # hide
+```
+
+![plot of the Rosenbrock function](obj_fun.png)
+
+### Explanation
+Every optimization problem consists of a set of _optimization variables_. In this case we create two variables. Additionally we assign _box constraints_ for each of them. In the next step we create two parameters for the problem with `@parameters`. While it is not needed to do this it makes it easier to `remake` the problem later with different values for the parameters. The _objective function_ is specified as well and finally everything is used to construct an `OptimizationSystem`.
 
+## Building and Solving the Optimization Problem
+Next the actual `OptimizationProblem` can be created. At this stage an initial guess `u0` for the optimization variables needs to be provided via map, using the symbols from before. Concrete values for the parameters of the system can also be provided or changed. However, if the parameters have default values assigned they are used automatically.
+```@example rosenbrock_2d
 u0 = [
-    x=>1.0
-    y=>2.0
+    x => 1.0
+    y => 2.0
 ]
 p = [
-    a => 6.0
-    b => 7.0
+    a => 1.0
+    b => 100.0
+]
+
+prob = OptimizationProblem(sys, u0, p, grad=true, hess=true)
+solve(prob, GradientDescent())
+```
+
+## Rosenbrock Function with Constraints
+```@example rosenbrock_2d_cstr
+using ModelingToolkit, Optimization, OptimizationOptimJL
+
+@variables begin
+    x, [bounds = (-2.0, 2.0)]
+    y, [bounds = (-1.0, 3.0)]
+end
+@parameters a = 1 b = 100
+loss = (a - x)^2 + b * (y - x^2)^2
+cons = [
+    x^2 + y^2 ≲ 1,
 ]
+@named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
+u0 = [
+    x => 1.0
+    y => 2.0
+]
+prob = OptimizationProblem(sys, u0, grad=true, hess=true)
+solve(prob, IPNewton())
+```
 
-prob = OptimizationProblem(sys,u0,p,grad=true,hess=true)
-solve(prob,Newton())
+A visualization of the objective function and the inequality constraint is depicted below.
+```@eval
+using Plots
+x = -2:0.01:2
+y = -1:0.01:3
+contour(x, y, (x,y) -> (1 - x)^2 + 100 * (y - x^2)^2, fill=true, color=:viridis, ratio=:equal, xlims=(-2, 2))
+contour!(x, y, (x, y) -> x^2 + y^2, levels=[1], color=:lightblue, line=4)
+savefig("obj_fun_c.png"); nothing # hide
 ```
 
+![plot of the Rosenbrock function with constraint](obj_fun_c.png)
+
+### Explanation
+Equality and inequality constraints can be added to the `OptimizationSystem`. An equality constraint can be specified via and `Equation`, e.g., `x^2 + y^2 ~ 1`. While inequality constraints via an `Inequality`, e.g., `x^2 + y^2 ≲ 1`. The syntax is here `\lesssim` and `\gtrsim`.
+
+## Nested Systems
 Needs more text but it's super cool and auto-parallelizes and sparsifies too.
 Plus you can hierarchically nest systems to have it generate huge
 optimization problems.

src/ModelingToolkit.jl

@@ -67,7 +67,15 @@ import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
 for fun in [:toexpr]
     @eval begin
         function $fun(eq::Equation; kw...)
-            Expr(:(=), $fun(eq.lhs; kw...), $fun(eq.rhs; kw...))
+            Expr(:call, :(==), $fun(eq.lhs; kw...), $fun(eq.rhs; kw...))
+        end
+
+        function $fun(ineq::Inequality; kw...)
+            if ineq.relational_op == Symbolics.leq
+                Expr(:call, :(<=), $fun(ineq.lhs; kw...), $fun(ineq.rhs; kw...))
+            else
+                Expr(:call, :(>=), $fun(ineq.lhs; kw...), $fun(ineq.rhs; kw...))
+            end
         end
 
         $fun(eqs::AbstractArray; kw...) = map(eq -> $fun(eq; kw...), eqs)
@@ -86,6 +94,7 @@ abstract type AbstractTimeDependentSystem <: AbstractSystem end
 abstract type AbstractTimeIndependentSystem <: AbstractSystem end
 abstract type AbstractODESystem <: AbstractTimeDependentSystem end
 abstract type AbstractMultivariateSystem <: AbstractSystem end
+abstract type AbstractOptimizationSystem <: AbstractTimeIndependentSystem end
 
 """
 $(TYPEDSIGNATURES)
@@ -139,6 +148,7 @@ include("systems/jumps/jumpsystem.jl")
 include("systems/nonlinear/nonlinearsystem.jl")
 include("systems/nonlinear/modelingtoolkitize.jl")
 
+include("systems/optimization/constraints_system.jl")
 include("systems/optimization/optimizationsystem.jl")
 
 include("systems/pde/pdesystem.jl")
@@ -179,7 +189,7 @@ export OptimizationProblem, OptimizationProblemExpr, constraints
 export AutoModelingToolkit
 export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem, DiscreteProblem
-export NonlinearSystem, OptimizationSystem
+export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
 export connect, @connector, Connection, Flow, Stream, instream
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist,
@@ -192,7 +202,7 @@ export Equation, ConstrainedEquation
 export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variables, independent_variable, states, parameters, equations, controls,
-       observed, structure, full_equations
+       observed, structure, full_equations, objective
 export structural_simplify, expand_connections, linearize, linearization_function
 export DiscreteSystem, DiscreteProblem
 

src/constants.jl

@@ -32,3 +32,15 @@ macro constants(xs...)
                           xs,
                           toconstant) |> esc
 end
+
+"""
+Substitute all `@constants` in the given expression
+"""
+function subs_constants(eqs)
+    consts = collect_constants(eqs)
+    if !isempty(consts)
+        csubs = Dict(c => getdefault(c) for c in consts)
+        eqs = substitute(eqs, csubs)
+    end
+    return eqs
+end

src/inputoutput.jl

@@ -343,14 +343,16 @@ function get_disturbance_system(dist::DisturbanceModel{<:ODESystem})
 end
 
 """
-    (f_oop, f_ip), augmented_sys, dvs, p = add_input_disturbance(sys, dist::DisturbanceModel)
+    (f_oop, f_ip), augmented_sys, dvs, p = add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing)
 
 Add a model of an unmeasured disturbance to `sys`. The disturbance model is an instance of [`DisturbanceModel`](@ref).
 
 The generated dynamics functions `(f_oop, f_ip)` will preserve any state and dynamics associated with disturbance inputs, but the disturbance inputs themselves will not be included as inputs to the generated function. The use case for this is to generate dynamics for state observers that estimate the influence of unmeasured disturbances, and thus require state variables for the disturbance model, but without disturbance inputs since the disturbances are not available for measurement.
 
 `dvs` will be the states of the simplified augmented system, consisting of the states of `sys` as well as the states of the disturbance model.
 
+For MIMO systems, all inputs to the system has to be specified in the argument `inputs`
+
 # Example
 The example below builds a double-mass model and adds an integrating disturbance to the input
 ```julia
@@ -391,18 +393,29 @@ dist = ModelingToolkit.DisturbanceModel(model.torque.tau.u, dmodel)
 ```
 `f_oop` will have an extra state corresponding to the integrator in the disturbance model. This state will not be affected by any input, but will affect the dynamics from where it enters, in this case it will affect additively from `model.torque.tau.u`.
 """
-function add_input_disturbance(sys, dist::DisturbanceModel)
+function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing)
     t = get_iv(sys)
     @variables d(t)=0 [disturbance = true]
-    @variables u(t)=0 [input = true]
+    @variables u(t)=0 [input = true] # New system input
     dsys = get_disturbance_system(dist)
 
+    if inputs === nothing
+        all_inputs = [u]
+    else
+        i = findfirst(isequal(dist.input), inputs)
+        if i === nothing
+            throw(ArgumentError("Input $(dist.input) indicated in the disturbance model was not found among inputs specified to add_input_disturbance"))
+        end
+        all_inputs = copy(inputs)
+        all_inputs[i] = u # The input where the disturbance acts is no longer an input, the new input is u
+    end
+
     eqs = [dsys.input.u[1] ~ d
            dist.input ~ u + dsys.output.u[1]]
 
     augmented_sys = ODESystem(eqs, t, systems = [sys, dsys], name = gensym(:outer))
 
-    (f_oop, f_ip), dvs, p = generate_control_function(augmented_sys, [u],
+    (f_oop, f_ip), dvs, p = generate_control_function(augmented_sys, all_inputs,
                                                       [d])
     (f_oop, f_ip), augmented_sys, dvs, p
 end

src/systems/diffeqs/abstractodesystem.jl

@@ -562,7 +562,7 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                           syms = $(Symbol.(states(sys))),
                           indepsym = $(QuoteNode(Symbol(get_iv(sys)))),
                           paramsyms = $(Symbol.(parameters(sys))),
-                          sparsity = $(jacobian_sparsity(sys)))
+                          sparsity = $sparsity ? $(jacobian_sparsity(sys)) : $nothing)
     end
     !linenumbers ? striplines(ex) : ex
 end