-
-
Notifications
You must be signed in to change notification settings - Fork 199
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
6 changed files
with
570 additions
and
178 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,26 +1,85 @@ | ||
| # 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)) | ||
| ``` | ||
|
|
||
| ### 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) | ||
| ``` | ||
|
|
||
| ### 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. |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,211 @@ | ||
| """ | ||
| $(TYPEDEF) | ||
| A constraint system of equations. | ||
| # Fields | ||
| $(FIELDS) | ||
| # Examples | ||
| ```julia | ||
| @variables x y z | ||
| @parameters a b c | ||
| cstr = [0 ~ a*(y-x), | ||
| 0 ~ x*(b-z)-y, | ||
| 0 ~ x*y - c*z | ||
| x^2 + y^2 ≲ 1] | ||
| @named ns = ConstraintsSystem(cstr, [x,y,z],[a,b,c]) | ||
| ``` | ||
| """ | ||
| struct ConstraintsSystem <: AbstractTimeIndependentSystem | ||
| """ | ||
| tag: a tag for the system. If two system have the same tag, then they are | ||
| structurally identical. | ||
| """ | ||
| tag::UInt | ||
| """Vector of equations defining the system.""" | ||
| constraints::Vector{Union{Equation, Inequality}} | ||
| """Unknown variables.""" | ||
| states::Vector | ||
| """Parameters.""" | ||
| ps::Vector | ||
| """Array variables.""" | ||
| var_to_name::Any | ||
| """Observed variables.""" | ||
| observed::Vector{Equation} | ||
| """ | ||
| Jacobian matrix. Note: this field will not be defined until | ||
| [`calculate_jacobian`](@ref) is called on the system. | ||
| """ | ||
| jac::RefValue{Any} | ||
| """ | ||
| Name: the name of the system. These are required to have unique names. | ||
| """ | ||
| name::Symbol | ||
| """ | ||
| systems: The internal systems | ||
| """ | ||
| systems::Vector{ConstraintsSystem} | ||
| """ | ||
| defaults: The default values to use when initial conditions and/or | ||
| parameters are not supplied in `ODEProblem`. | ||
| """ | ||
| defaults::Dict | ||
| """ | ||
| type: type of the system | ||
| """ | ||
| connector_type::Any | ||
| """ | ||
| metadata: metadata for the system, to be used by downstream packages. | ||
| """ | ||
| metadata::Any | ||
| """ | ||
| tearing_state: cache for intermediate tearing state | ||
| """ | ||
| tearing_state::Any | ||
| """ | ||
| substitutions: substitutions generated by tearing. | ||
| """ | ||
| substitutions::Any | ||
|
|
||
| function ConstraintsSystem(tag, constraints, states, ps, var_to_name, observed, jac, | ||
| name, | ||
| systems, | ||
| defaults, connector_type, metadata = nothing, | ||
| tearing_state = nothing, substitutions = nothing; | ||
| checks::Union{Bool, Int} = true) | ||
| if checks == true || (checks & CheckUnits) > 0 | ||
| all_dimensionless([states; ps]) || check_units(constraints) | ||
| end | ||
| new(tag, constraints, states, ps, var_to_name, observed, jac, name, systems, | ||
| defaults, | ||
| connector_type, metadata, tearing_state, substitutions) | ||
| end | ||
| end | ||
|
|
||
| equations(sys::ConstraintsSystem) = constraints(sys) # needed for Base.show | ||
|
|
||
| function ConstraintsSystem(constraints, states, ps; | ||
| observed = [], | ||
| name = nothing, | ||
| default_u0 = Dict(), | ||
| default_p = Dict(), | ||
| defaults = _merge(Dict(default_u0), Dict(default_p)), | ||
| systems = ConstraintsSystem[], | ||
| connector_type = nothing, | ||
| continuous_events = nothing, # this argument is only required for ODESystems, but is added here for the constructor to accept it without error | ||
| discrete_events = nothing, # this argument is only required for ODESystems, but is added here for the constructor to accept it without error | ||
| checks = true, | ||
| metadata = nothing) | ||
| continuous_events === nothing || isempty(continuous_events) || | ||
| throw(ArgumentError("ConstraintsSystem does not accept `continuous_events`, you provided $continuous_events")) | ||
| discrete_events === nothing || isempty(discrete_events) || | ||
| throw(ArgumentError("ConstraintsSystem does not accept `discrete_events`, you provided $discrete_events")) | ||
|
|
||
| name === nothing && | ||
| throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro")) | ||
|
|
||
| cstr = value.(Symbolics.canonical_form.(scalarize(constraints))) | ||
| states′ = value.(scalarize(states)) | ||
| ps′ = value.(scalarize(ps)) | ||
|
|
||
| if !(isempty(default_u0) && isempty(default_p)) | ||
| Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", | ||
| :ConstraintsSystem, force = true) | ||
| end | ||
| sysnames = nameof.(systems) | ||
| if length(unique(sysnames)) != length(sysnames) | ||
| throw(ArgumentError("System names must be unique.")) | ||
| end | ||
|
|
||
| jac = RefValue{Any}(EMPTY_JAC) | ||
| defaults = todict(defaults) | ||
| defaults = Dict(value(k) => value(v) for (k, v) in pairs(defaults)) | ||
|
|
||
| var_to_name = Dict() | ||
| process_variables!(var_to_name, defaults, states′) | ||
| process_variables!(var_to_name, defaults, ps′) | ||
| isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed)) | ||
|
|
||
| ConstraintsSystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), | ||
| cstr, states, ps, var_to_name, observed, jac, name, systems, | ||
| defaults, | ||
| connector_type, metadata, checks = checks) | ||
| end | ||
|
|
||
| function calculate_jacobian(sys::ConstraintsSystem; sparse = false, simplify = false) | ||
| cache = get_jac(sys)[] | ||
| if cache isa Tuple && cache[2] == (sparse, simplify) | ||
| return cache[1] | ||
| end | ||
|
|
||
| lhss = generate_canonical_form_lhss(sys) | ||
| vals = [dv for dv in states(sys)] | ||
| if sparse | ||
| jac = sparsejacobian(lhss, vals, simplify = simplify) | ||
| else | ||
| jac = jacobian(lhss, vals, simplify = simplify) | ||
| end | ||
| get_jac(sys)[] = jac, (sparse, simplify) | ||
| return jac | ||
| end | ||
|
|
||
| function generate_jacobian(sys::ConstraintsSystem, vs = states(sys), ps = parameters(sys); | ||
| sparse = false, simplify = false, kwargs...) | ||
| jac = calculate_jacobian(sys, sparse = sparse, simplify = simplify) | ||
| return build_function(jac, vs, ps; kwargs...) | ||
| end | ||
|
|
||
| function calculate_hessian(sys::ConstraintsSystem; sparse = false, simplify = false) | ||
| lhss = generate_canonical_form_lhss(sys) | ||
| vals = [dv for dv in states(sys)] | ||
| if sparse | ||
| hess = [sparsehessian(lhs, vals, simplify = simplify) for lhs in lhss] | ||
| else | ||
| hess = [hessian(lhs, vals, simplify = simplify) for lhs in lhss] | ||
| end | ||
| return hess | ||
| end | ||
|
|
||
| function generate_hessian(sys::ConstraintsSystem, vs = states(sys), ps = parameters(sys); | ||
| sparse = false, simplify = false, kwargs...) | ||
| hess = calculate_hessian(sys, sparse = sparse, simplify = simplify) | ||
| return build_function(hess, vs, ps; kwargs...) | ||
| end | ||
|
|
||
| function generate_function(sys::ConstraintsSystem, dvs = states(sys), ps = parameters(sys); | ||
| kwargs...) | ||
| lhss = generate_canonical_form_lhss(sys) | ||
| pre, sol_states = get_substitutions_and_solved_states(sys) | ||
|
|
||
| func = build_function(lhss, value.(dvs), value.(ps); postprocess_fbody = pre, | ||
| states = sol_states, kwargs...) | ||
|
|
||
| cstr = constraints(sys) | ||
| lcons = fill(-Inf, length(cstr)) | ||
| ucons = zeros(length(cstr)) | ||
| lcons[findall(Base.Fix2(isa, Equation), cstr)] .= 0.0 | ||
|
|
||
| return func, lcons, ucons | ||
| end | ||
|
|
||
| function jacobian_sparsity(sys::ConstraintsSystem) | ||
| lhss = generate_canonical_form_lhss(sys) | ||
| jacobian_sparsity(lhss, states(sys)) | ||
| end | ||
|
|
||
| function hessian_sparsity(sys::ConstraintsSystem) | ||
| lhss = generate_canonical_form_lhss(sys) | ||
| [hessian_sparsity(eq, states(sys)) for eq in lhss] | ||
| end | ||
|
|
||
| """ | ||
| Convert the system of equalities and inequalities into a canonical form: | ||
| h(x) = 0 | ||
| g(x) <= 0 | ||
| """ | ||
| function generate_canonical_form_lhss(sys) | ||
| lhss = [Symbolics.canonical_form(eq).lhs for eq in constraints(sys)] | ||
| end |
Oops, something went wrong.