Implement constant_over_collocation

docs/src/examples/Optimal Control/pandemic_control.jl

@@ -13,11 +13,10 @@
 # This model is formalized as:
 # ```math 
 # \begin{gathered}
-# \frac{ds(t)}{dt} = (u(t) - 1)\beta si(t) \\
-# \frac{de(t)}{dt} = (1 - u(t))\beta si(t) - \xi e(t) \\
+# \frac{ds(t)}{dt} = (u(t) - 1)\beta s(t)i(t) \\
+# \frac{de(t)}{dt} = (1 - u(t))\beta s(t)i(t) - \xi e(t) \\
 # \frac{di(t)}{dt} = \xi e(t) - \gamma i(t)\\
 # \frac{dr(t)}{dt} = \gamma i(t) \\
-# si(t) = s(t) i(t)
 # \end{gathered}
 # ```
 # where ``s(t)`` is the susceptible population, ``e(t)`` is the exposed population, 
@@ -36,11 +35,10 @@
 # ```math 
 # \begin{aligned}
 # &&\min_{} &&& \int_{t \in \mathcal{D}_{t}} u(t) dt \\
-# && \text{s.t.} &&& \frac{\partial s(t, \xi)}{\partial t} = (u(t) - 1)\beta si(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
-# &&&&& \frac{\partial e(t, \xi)}{\partial t} = (1 - u(t))\beta si(t, \xi) - \xi e(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
+# && \text{s.t.} &&& \frac{\partial s(t, \xi)}{\partial t} = (u(t) - 1)\beta s(t, \xi)i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
+# &&&&& \frac{\partial e(t, \xi)}{\partial t} = (1 - u(t))\beta s(t, \xi)i(t, \xi) - \xi e(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
 # &&&&& \frac{\partial i(t, \xi)}{\partial t} = \xi e(t, \xi) - \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
 # &&&&& \frac{\partial r(t, \xi)}{\partial t} = \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
-# &&&&& si(t, \xi) = s(t, \xi) i(t, \xi), && \forall \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
 # &&&&& s(0, \xi) = s_0, e(0, \xi) = e_0, i(0, \xi) = i_0, r(0, \xi) = r_0, && \forall \xi \in \mathcal{D}_{\xi} \\
 # &&&&& i(t, \xi) \leq i_{max}, && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
 # &&&&& u(t) \in [0, 0.8] \\
@@ -100,7 +98,7 @@ model = InfiniteModel(Ipopt.Optimizer);
 # - specify that the number of random scenarios should equal `num_samples`
 # - add `extra_ts` as extra time points
 @infinite_parameter(model, t ∈ [t0, tf], num_supports = 51, 
-                    derivative_method = OrthogonalCollocation(2))
+                    derivative_method = OrthogonalCollocation(3))
 @infinite_parameter(model, ξ ~ Uniform(ξ_min, ξ_max), num_supports = num_samples)
 add_supports(t, extra_ts)
 
@@ -110,13 +108,11 @@ add_supports(t, extra_ts)
 # - ``e(t, \xi) \geq 0``
 # - ``i(t, \xi) \geq 0``
 # - ``r(t, \xi) \geq 0``
-# - ``si(t, \xi)``
 # - ``0 \leq u(t) \leq 0.8``
 @variable(model, s ≥ 0, Infinite(t, ξ))
 @variable(model, e ≥ 0, Infinite(t, ξ))
 @variable(model, i ≥ 0, Infinite(t, ξ))
 @variable(model, r ≥ 0, Infinite(t, ξ))
-@variable(model, si, Infinite(t, ξ))
 @variable(model, 0 ≤ u ≤ 0.8, Infinite(t), start = 0.2)
 
 # ## Objective Definition
@@ -137,7 +133,6 @@ add_supports(t, extra_ts)
 # &&& \frac{\partial e(t, \xi)}{\partial t} = (1 - u(t))\beta si(t, \xi) - \xi e(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
 # &&& \frac{\partial i(t, \xi)}{\partial t} = \xi e(t, \xi) - \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
 # &&& \frac{\partial r(t, \xi)}{\partial t} = \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
-# &&& si(t, \xi) = s(t, \xi) i(t, \xi), && \forall \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi}, \\
 # \end{aligned}
 # ```
 # and the infection limit constraint:
@@ -151,15 +146,21 @@ add_supports(t, extra_ts)
 @constraint(model, r(0, ξ) == r0)
 
 ## Define the SEIR equations
-@constraint(model, s_constr, ∂(s, t) == -(1 - u) * β * si)
-@constraint(model, e_constr, ∂(e, t) == (1 - u) * β * si - ξ * e)
+@constraint(model, s_constr, ∂(s, t) == -(1 - u) * β * s * i)
+@constraint(model, e_constr, ∂(e, t) == (1 - u) * β * s * i - ξ * e)
 @constraint(model, i_constr, ∂(i, t) == ξ * e - γ * i)
 @constraint(model, r_constr, ∂(r, t) == γ * i)
-@constraint(model, si == s * i)
 
 ## Define the infection rate limit
 @constraint(model, imax_constr, i ≤ i_max)
 
+# ## Adjust Degrees of Freedom
+# Since we are using [`OrthogonalCollocation`](@ref), this introduces additional collocation points 
+# for `t`. These in turn, will artificially increase the degrees of freedom for `u`. Thus, 
+# we will treat `u` as a piecewise constant function that is held constant over the internal 
+# collocation nodes via [`constant_over_collocation`](@ref constant_over_collocation(::InfiniteVariableRef, ::GeneralVariableRef)):
+constant_over_collocation(u, t)
+
 # ## Display the Infinite Model
 # Let's display `model` now that it is fully defined:
 print(model)

docs/src/guide/derivative.md

@@ -37,35 +37,35 @@ over finite elements using 3 nodes. We'll come back to this just a little furthe
 below to more fully describe the various methods we can use. 
 
 First, let's discuss how to define derivatives in `InfiniteOpt.jl`. Principally, 
-this is accomplished via [`@deriv`](@ref) which will operate on a particular 
+this is accomplished via [`deriv`](@ref) which will operate on a particular 
 `InfiniteOpt` expression (containing parameters, variables, and/or measures) with 
 respect to infinite parameters specified with their associated orders. Behind the 
 scenes all the appropriate calculus will be applied, creating derivative variables 
 as needed. For example, we can define the following:
 ```jldoctest deriv_basic 
-julia> d1 = @deriv(y, t)
+julia> d1 = deriv(y, t)
 ∂/∂t[y(t, ξ)]
 
-julia> d2 = @deriv(y, t, ξ)
+julia> d2 = ∂(y, t, ξ)
 ∂/∂ξ[∂/∂t[y(t, ξ)]]
 
-julia> d3 = @∂(q, t^2)
+julia> d3 = @∂(q, t^2) # the macro version allows the `t^2` syntax
 ∂/∂t[∂/∂t[q(t)]]
 
-julia> d_expr = @deriv(y * q - 2t, t)
+julia> d_expr = deriv(y * q - 2t, t)
 ∂/∂t[y(t, ξ)]*q(t) + y(t, ξ)*∂/∂t[q(t)] - 2
 ```
 Thus, we can define derivatives in a variety of forms according to the problem at 
 hand. The last example even shows how the product rule is correctly applied. 
 
 !!! note 
-    For convenience in making more compact code we provide [`∂`](@ref) and 
-    [`@∂`](@ref) as wrappers for [`deriv`](@ref) and [`@deriv`](@ref), respectively.
+    For convenience in making more compact code we provide [`∂`](@ref) 
+    as a wrapper for [`deriv`](@ref).
 
 Also, notice that the appropriate analytic calculus is applied to infinite 
 parameters. For example, we could also compute:
 ```jldoctest deriv_basic 
-julia> @deriv(3t^2 - 2t, t)
+julia> deriv(3t^2 - 2t, t)
 6 t - 2
 ```
 
@@ -109,6 +109,26 @@ julia> set_all_derivative_methods(model, FiniteDifference(Forward()))
 
 ```
 
+!!! note
+    When [`OrthogonalCollocation`](@ref) is used, additional degrees of freedom can be 
+    artificially introduced to infinite variables that share the same infinite parameter. 
+    For instance, this occurs with control variables in optimal control problems. To address 
+    this, [`constant_over_collocation`](@ref constant_over_collocation(::InfiniteVariableRef, ::GeneralVariableRef)) 
+    should be called on the appropriate variables. For example:
+    ```jldoctest collocation; setup = :(using InfiniteOpt; model = InfiniteModel()) 
+    @infinite_parameter(model, t in [0, 1], derivative_method = OrthogonalCollocation(3))
+    @variable(model, y_state, Infinite(t))
+    @variable(model, y_control, Infinite(t))
+    @constraint(model, ∂(y_state, t) == y_state^2)
+    @constraint(model, y_state(0) == 0)
+    constant_over_collocation(y_control, t)
+
+    # output
+
+    ```
+    where we use `constant_over_collocation` to hold `y_control` constant over each finite 
+    element (i.e., constant for each internal collocation point). 
+
 !!! warning
     `InfiniteOpt` does not ensure proper boundary conditions are provided by the 
     user. Thus, it is imperative that the user ensure these are provided appropriately 
@@ -388,6 +408,11 @@ orthogonal collocation over finite elements, this
 [page](http://apmonitor.com/do/index.php/Main/OrthogonalCollocation) provides a 
 good reference.
 
+The addition of internal collocation supports by `OrthogonalCollocation` will increase 
+the degrees of freedom for infinite variables that are not used by derivatives (e.g., 
+control variables). To prevent this, we use [`constant_over_collocation`](@ref) on any 
+such infinite variables to hold them constant over internal collocation nodes. 
+
 Other methods can be employed via user-defined extensions. Please visit our 
 Extensions page for more information.
 

docs/src/guide/measure.md

@@ -190,7 +190,7 @@ The arguments of [`DiscreteMeasureData`](@ref) are parameter, coefficients, and
 supports. The default weight function is ``w(\tau) = 1`` for 
 any ``\tau``, which can be overwritten by the keyword argument `weight_function`. 
 The `weight_function` should take a function that returns a number for any 
-value that is well defined for the integrated infinite parameter. The data type 
+value that is well-defined for the integrated infinite parameter. The data type 
 is [`DiscreteMeasureData`](@ref), which is a subtype of the abstract data type 
 [`AbstractMeasureData`](@ref).
 

docs/src/manual/expression.md

@@ -154,6 +154,7 @@ JuMP.is_integer(::GeneralVariableRef)
 JuMP.set_integer(::GeneralVariableRef)
 JuMP.IntegerRef(::GeneralVariableRef)
 JuMP.unset_integer(::GeneralVariableRef)
+constant_over_collocation(::GeneralVariableRef, ::GeneralVariableRef)
 ```
 
 ## Developer Internal Methods

docs/src/manual/transcribe.md

@@ -17,6 +17,7 @@ InfiniteOpt.TranscriptionOpt.transcribe_measures!
 InfiniteOpt.TranscriptionOpt.transcribe_objective!
 InfiniteOpt.TranscriptionOpt.transcribe_constraints!
 InfiniteOpt.TranscriptionOpt.transcribe_derivative_evaluations!
+InfiniteOpt.TranscriptionOpt.transcribe_variable_collocation_restictions!
 InfiniteOpt.TranscriptionOpt.build_transcription_model!
 InfiniteOpt.add_point_variable(::JuMP.Model,::InfiniteOpt.GeneralVariableRef,::Vector{Float64},::Val{:TransData})
 InfiniteOpt.add_semi_infinite_variable(::JuMP.Model,::InfiniteOpt.SemiInfiniteVariable,::Val{:TransData})

docs/src/manual/variable.md

@@ -144,6 +144,7 @@ JuMP.delete(::InfiniteModel, ::DecisionVariableRef)
 ```@docs
 set_start_value_function(::InfiniteVariableRef, ::Union{Real, Function})
 reset_start_value_function(::InfiniteVariableRef)
+constant_over_collocation(::InfiniteVariableRef, ::GeneralVariableRef)
 ```
 
 ### Semi-Infinite

docs/src/tutorials/quick_start.md

@@ -95,7 +95,7 @@ this problem we have the time domain ``t \in \mathcal{D}_t`` and the random doma
 ``\xi \in \mathcal{D}_\xi`` where ``\xi \sim \mathcal{N}(\mu, \sigma^2)``:
 ```jldoctest quick
 julia> @infinite_parameter(model, t in [0, 60], num_supports = 61, 
-                           derivative_method = OrthogonalCollocation(2))
+                           derivative_method = OrthogonalCollocation(3))
 t
 
 julia> @infinite_parameter(model, ξ ~ Normal(μ, σ^2), num_supports = 10)
@@ -147,7 +147,7 @@ And data, a 4-element Vector{GeneralVariableRef}:
  y[3](ξ)
  y[4](ξ)
 ```
-Notice that we specifying the initial guess for all of them via `start`. We also 
+Notice that we specify the initial guess for all of them via `start`. We also 
 can symbolically define variable conditions like the lower bound on `y`.
 
 That does it for this example, but other problems might also employ the following:
@@ -212,7 +212,7 @@ And data, a 2-element Vector{InfOptConstraintRef}:
  c2[2] : ξ*∂/∂t[v[2](t, ξ)] - u[2](t) = 0, ∀ t ∈ [0, 60], ξ ~ Normal
 ```
 
-Finally, we can define our last 2 constraints:
+Next, we can define our last 2 constraints:
  ```jldoctest quick
 julia> @constraint(model, c3[w in W], y[w] == sum((x[i](tw[w], ξ) - p[i, w])^2 for i in I))
 1-dimensional DenseAxisArray{InfOptConstraintRef,1,...} with index sets:
@@ -228,18 +228,26 @@ c4 : 𝔼{ξ}[y[1](ξ) + y[2](ξ) + y[3](ξ) + y[4](ξ)] - ϵ ≤ 0
 ```
 Notice we are able to invoke an expectation simply by calling [`expect`](@ref).
 
+Finally, to address any unwanted degrees of freedom introduced by internal collocation 
+nodes with [`OrthogonalCollocation`](@ref). We should call [`constant_over_collocation`](@ref constant_over_collocation(::InfiniteVariableRef, ::GeneralVariableRef)) 
+on any control variables:
+```jldoctest quick
+julia> constant_over_collocation.(u, t);
+
+``` 
+
 That's it, now we have our problem defined in `InfiniteOpt`!
 
 ## Solution & Queries
 ### Optimize 
 Now that our model is defined, let's optimize it via [`optimize!`](@ref):
-```julia-repl
+```jldoctest quick; setup = :(set_optimizer_attribute(model, "print_level", 0))
 julia> optimize!(model)
 
 ```
 We can check the solution status via 
 [`termination_status`](@ref JuMP.termination_status(::InfiniteModel)):
-```jldoctest quick; setup = :(set_optimizer_attribute(model, "print_level", 0); optimize!(model))
+```jldoctest quick
 julia> termination_status(model)
 LOCALLY_SOLVED::TerminationStatusCode = 4
 ```
@@ -287,7 +295,7 @@ model = InfiniteModel(Ipopt.Optimizer)
 # INITIALIZE THE PARAMETERS
 @finite_parameter(model, ϵ == 10)
 @infinite_parameter(model, t in [0, 60], num_supports = 61, 
-                    derivative_method = OrthogonalCollocation(2))
+                    derivative_method = OrthogonalCollocation(3))
 @infinite_parameter(model, ξ ~ Normal(μ, σ^2), num_supports = 10)
 
 # INITIALIZE THE VARIABLES
@@ -309,6 +317,9 @@ model = InfiniteModel(Ipopt.Optimizer)
 @constraint(model, c3[w in W], y[w] == sum((x[i](tw[w], ξ) - p[i, w])^2 for i in I))
 @constraint(model, c4, expect(sum(y[w] for w in W), ξ) <= ϵ)
 
+# ADJUST DEGREES OF FREEDOM FOR CONTROL VARIABLES
+constant_over_collocation.(u, t)
+
 # SOLVE THE MODEL
 optimize!(model)
 

src/TranscriptionOpt/model.jl

@@ -907,7 +907,7 @@ function support_index_iterator(
     raw_supps = parameter_supports(model)
     lens = map(i -> length(i), raw_supps)
     # prepare the indices of each support combo
-    # note that the actual supports are afrom 1:length-1 and the placeholders are at the ends
+    # note that the actual supports are from 1:length-1 and the placeholders are at the ends
     return CartesianIndices(ntuple(i -> i in obj_nums ? (1:lens[i]-1) : (lens[i]:lens[i]),
                                    length(raw_supps)))
 end

src/TranscriptionOpt/transcribe.jl

@@ -782,7 +782,7 @@ function transcribe_constraints!(
 end
 
 ################################################################################
-#                      DERIVATIVE CONSTRAINT TRANSCRIPTION METHODS
+#                    DERIVATIVE CONSTRAINT TRANSCRIPTION METHODS
 ################################################################################
 """
     transcribe_derivative_evaluations!(trans_model::JuMP.Model, 
@@ -791,7 +791,7 @@ end
 Generate the auxiliary derivative evaluation equations and transcribe them 
 appropriately for all the derivatives in `inf_model`. These are in turn added to 
 `trans_model`. Note that no mapping information is recorded since the InfiniteModel 
-won't have any constraints that correspond to these equations. Also Note that
+won't have any constraints that correspond to these equations. Also note that
 the variables and measures must all first be transcripted (e.g., via
 [`transcribe_derivative_variables!`](@ref)).
 """
@@ -827,6 +827,52 @@ function transcribe_derivative_evaluations!(
     return
 end
 
+"""
+    transcribe_variable_collocation_restictions!(trans_model::JuMP.Model, 
+                                                 inf_model::InfiniteOpt.InfiniteModel)::Nothing
+
+Add constraints to `trans_model` that make infinite variables constant over collocation 
+points following the calls made to [`InfiniteOpt.constant_over_collocation`](@ref). Note that 
+[`set_parameter_supports`](@ref) and [`transcribe_infinite_variables!`](@ref) must be called first.
+"""
+function transcribe_variable_collocation_restictions!(
+    trans_model::JuMP.Model, 
+    inf_model::InfiniteOpt.InfiniteModel
+    )
+    data = transcription_data(trans_model)
+    set = MOI.EqualTo(0.0)
+    for (pidx, vidxs) in inf_model.piecewise_vars
+        pref = InfiniteOpt._make_variable_ref(inf_model, pidx)
+        if !InfiniteOpt.has_generative_supports(pref)
+            continue
+        end
+        obj_num = InfiniteOpt._object_number(pref)
+        supps = reverse!(data.supports[obj_num][1:end-1])
+        labels = reverse!(data.support_labels[obj_num][1:end-1])
+        @assert any(l -> l <: InfiniteOpt.PublicLabel, first(labels))
+        v_manip = GeneralVariableRef(inf_model, -1, IndependentParameterIndex) # placeholder
+        for vidx in vidxs
+            vref = InfiniteOpt._make_variable_ref(inf_model, vidx)
+            obj_nums = filter(!isequal(obj_num), InfiniteOpt._object_numbers(vref))
+            supp_indices = support_index_iterator(trans_model, obj_nums)
+            for (s, ls) in zip(supps, labels)
+                if any(l -> l <: InfiniteOpt.PublicLabel, ls)
+                    v_manip = InfiniteOpt.make_reduced_expr(vref, pref, s, trans_model)
+                else
+                    inf_expr = v_manip - InfiniteOpt.make_reduced_expr(vref, pref, s, trans_model)
+                    for i in supp_indices 
+                        raw_supp = index_to_support(trans_model, i)
+                        new_expr = transcription_expression(trans_model, inf_expr, raw_supp)
+                        trans_constr = JuMP.build_constraint(error, new_expr, set)
+                        JuMP.add_constraint(trans_model, trans_constr)
+                    end
+                end
+            end
+        end
+    end
+    return
+end
+
 ################################################################################
 #                      INFINITEMODEL TRANSCRIPTION METHODS
 ################################################################################
@@ -881,6 +927,8 @@ function build_transcription_model!(
     transcribe_constraints!(trans_model, inf_model)
     # define the derivative evaluation constraints
     transcribe_derivative_evaluations!(trans_model, inf_model)
+    # define constraints for variables that are constant over collocation points
+    transcribe_variable_collocation_restictions!(trans_model, inf_model)
     return
 end