Merge branch 'master' into func-affect

Author: dodoplus

Project.toml

@@ -30,7 +30,6 @@ NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
-Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -66,7 +65,6 @@ NaNMath = "0.3, 1"
 NonlinearSolve = "0.3.8"
 RecursiveArrayTools = "2.3"
 Reexport = "0.2, 1"
-Requires = "1.0"
 RuntimeGeneratedFunctions = "0.4.3, 0.5"
 SciMLBase = "1.26.2"
 Setfield = "0.7, 0.8, 1"
@@ -79,9 +77,13 @@ Unitful = "1.1"
 julia = "1.6"
 
 [extras]
+AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -93,4 +95,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "ForwardDiff", "Optimization", "OptimizationOptimJL", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]
+test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

docs/pages.jl

@@ -8,7 +8,6 @@ pages = [
                                          "tutorials/nonlinear.md",
                                          "tutorials/optimization.md",
                                          "tutorials/stochastic_diffeq.md",
-                                         "tutorials/nonlinear_optimal_control.md",
                                          "tutorials/parameter_identifiability.md"],
     "ModelingToolkitize Tutorials" => Any["mtkitize_tutorials/modelingtoolkitize.md",
                                           "mtkitize_tutorials/modelingtoolkitize_index_reduction.md",
@@ -25,7 +24,6 @@ pages = [
                           "systems/JumpSystem.md",
                           "systems/NonlinearSystem.md",
                           "systems/OptimizationSystem.md",
-                          "systems/ControlSystem.md",
                           "systems/PDESystem.md"],
     "comparison.md",
     "internals.md",

docs/src/systems/ControlSystem.md

@@ -60,8 +60,6 @@ import DiffEqBase: @add_kwonly
 
 import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
 
-using Requires
-
 for fun in [:toexpr]
     @eval begin
         function $fun(eq::Equation; kw...)
@@ -137,8 +135,6 @@ include("systems/nonlinear/modelingtoolkitize.jl")
 
 include("systems/optimization/optimizationsystem.jl")
 
-include("systems/control/controlsystem.jl")
-
 include("systems/pde/pdesystem.jl")
 
 include("systems/sparsematrixclil.jl")
@@ -173,13 +169,11 @@ export AutoModelingToolkit
 export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem
-export ControlSystem
 export alias_elimination, flatten
 export connect, @connector, Connection, Flow, Stream, instream
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist,
        tunable_parameters, isirreducible
 export ode_order_lowering, dae_order_lowering, liouville_transform
-export runge_kutta_discretize
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation

docs/src/tutorials/nonlinear_optimal_control.md

@@ -184,6 +184,53 @@ function complete(g::BipartiteGraph{I}) where {I}
     BipartiteGraph(g.ne, g.fadjlist, badjlist)
 end
 
+# Matrix whose only purpose is to pretty-print the bipartite graph
+struct BipartiteAdjacencyList
+    u::Union{Vector{Int}, Nothing}
+end
+function Base.show(io::IO, l::BipartiteAdjacencyList)
+    if l.u === nothing
+        printstyled(io, '⋅', color = :light_black)
+    elseif isempty(l.u)
+        printstyled(io, '∅', color = :light_black)
+    else
+        print(io, l.u)
+    end
+end
+
+struct Label
+    s::String
+end
+Base.show(io::IO, l::Label) = print(io, l.s)
+
+struct BipartiteGraphPrintMatrix <:
+       AbstractMatrix{Union{Label, Int, BipartiteAdjacencyList}}
+    bpg::BipartiteGraph
+end
+Base.size(bgpm::BipartiteGraphPrintMatrix) = (max(nsrcs(bgpm.bpg), ndsts(bgpm.bpg)) + 1, 3)
+function Base.getindex(bgpm::BipartiteGraphPrintMatrix, i::Integer, j::Integer)
+    checkbounds(bgpm, i, j)
+    if i == 1
+        return (Label.(("#", "src", "dst")))[j]
+    elseif j == 1
+        return i - 1
+    elseif j == 2
+        return BipartiteAdjacencyList(i - 1 <= nsrcs(bgpm.bpg) ?
+                                      𝑠neighbors(bgpm.bpg, i - 1) : nothing)
+    elseif j == 3
+        return BipartiteAdjacencyList(i - 1 <= ndsts(bgpm.bpg) ?
+                                      𝑑neighbors(bgpm.bpg, i - 1) : nothing)
+    else
+        @assert false
+    end
+end
+
+function Base.show(io::IO, b::BipartiteGraph)
+    print(io, "BipartiteGraph with (", length(b.fadjlist), ", ",
+          isa(b.badjlist, Int) ? b.badjlist : length(b.badjlist), ") (𝑠,𝑑)-vertices\n")
+    Base.print_matrix(io, BipartiteGraphPrintMatrix(b))
+end
+
 """
 ```julia
 Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}

src/ModelingToolkit.jl

@@ -37,17 +37,57 @@ end
 
 function tear_graph_modia(structure::SystemStructure; varfilter = v -> true,
                           eqfilter = eq -> true)
+    # It would be possible here to simply iterate over all variables and attempt to
+    # use tearEquations! to produce a matching that greedily selects the minimal
+    # number of torn variables. However, we can do this process faster if we first
+    # compute the strongly connected components. In the absence of cycles and
+    # non-solvability, a maximal matching on the original graph will give us an
+    # optimal assignment. However, even with cycles, we can use the maximal matching
+    # to give us a good starting point for a good matching and then proceed to
+    # reverse edges in each scc to improve the solution. Note that it is possible
+    # to have optimal solutions that cannot be found by this process. We will not
+    # find them here [TODO: It would be good to have an explicit example of this.]
+
     @unpack graph, solvable_graph = structure
     var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter))
     var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
 
+    # Here, we're using a maximal matching on the post-pantelides system to find
+    # the strongly connected components of the system (of variables that depend
+    # on each other). The strongly connected components are unique, however, the
+    # maximal matching itself is not. Every maximal matching gives rise to the
+    # same set of strongly connected components, but the associated equations need
+    # not be the same. In the absence of solvability constraints, this may be a
+    # small issue, but here it is possible that an equation got assigned to an
+    # scc that cannot actually use it for solving a variable, but still precludes
+    # another scc from using it. To avoid this, we delete any assignments that
+    # are not in the solvable graph and extend the set of considered eqauations
+    # below.
+    for var in ndsts(solvable_graph)
+        var_eq_matching[var] === unassigned && continue
+        if !(BipartiteEdge(var, var_eq_matching[var]) in solvable_graph)
+            var_eq_matching[var] = unassigned
+        end
+    end
+
     for vars in var_sccs
         filtered_vars = filter(varfilter, vars)
         ieqs = Int[var_eq_matching[v]
                    for v in filtered_vars if var_eq_matching[v] !== unassigned]
         for var in vars
             var_eq_matching[var] = unassigned
         end
+        for var in filtered_vars
+            # Add any equations that we may not have been able to use earlier to see
+            # if a different matching may have been possible.
+            for eq′ in 𝑑neighbors(solvable_graph, var)
+                eqfilter(eq′) || continue
+                eq′ in ieqs && continue
+                if invview(var_eq_matching)[eq′] === unassigned
+                    push!(ieqs, eq′)
+                end
+            end
+        end
         tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars)
     end