Edit on GitHub

Quick Start Guide

Below we exemplify and briefly explain the very basics behind defining and solving an infinite-dimensional optimization problem in InfiniteOpt. Please refer to the Guide on our subsequent pages for more complete information. The Basic Usage sections on the each guide page are good places to start from. Also, the syntax of InfiniteOpt is inspired by JuMP thus we recommend new users that haven't used JuMP, first consult their tutorials starting here.

Preliminaries

Software Setup

First, we need to make sure everything is installed. This will include:

  • installing Julia
  • installing the InfiniteOpt.jl, JuMP.jl, and Distributions.jl packages
  • installing wanted optimizers e.g., Ipopt.jl and HiGHS.jl

See Installation for more information.

Problem Formulation

Now we need to formulate the problem we want to solve mathematically. For example, let's define a simple optimal control model:

\[\begin{aligned} &&\underset{x_i(t, \xi), v_i(t, \xi), y_w(\xi), u_i(t)}{\text{min}} &&& \int_{t \in \mathcal{D}_t} \sum_{i \in I} u_i^2(t) dt \\ &&\text{s.t.} &&& x_i(0, \xi) = x0_i, && \forall i \in I, \xi \in \mathcal{D}_\xi\\ &&&&& v_i(0, \xi) = v0_i, && \forall i \in I, \xi \in \mathcal{D}_\xi \\ &&&&& \frac{\partial x_i(t, \xi)}{\partial t} = v_i(t, \xi), && \forall i \in I, t \in \mathcal{D}_t, \xi \in \mathcal{D}_\xi\\ &&&&& \xi\frac{\partial v_i(t, \xi)}{\partial t} = u_i(t), && \forall i \in I, t \in \mathcal{D}_t, \xi \in \mathcal{D}_\xi\\ &&&&& y_{w}(\xi) = \sum_{i \in I}(x_i(t_w, \xi) - p_{iw})^2, && \forall w \in W, \xi \in \mathcal{D}_\xi \\ &&&&& y_{w}(\xi) \geq 0, && \forall w \in W, \xi \in \mathcal{D}_\xi \\ &&&&& \mathbb{E}_{\xi}\left[\sum_{w \in W} y_w(\xi) \right] \leq \epsilon \\ &&&&& \xi \sim \mathcal{N}(\mu, \sigma^2) \\ &&&&& t \in \mathcal{D}_t \end{aligned}\]

Notice this model is both dynamic with time $t$ and random with respect to $\xi$.

Parameter Specification

Before moving on we'll need to define the necessary constants and problem parameters. Thus, continuing with our example we define the following in our Julia session (these could also be put into a script as is shown at the bottom of this page):

julia> μ = 1; σ = 0.2; # set the distribution parameters 

julia> x0 = [0, 0]; v0 = [0, 0]; # set the initial conditions

julia> p = [1 4 6 1; 1 3 0 1]; tw = [0, 25, 50, 60]; # set waypoint specifications

julia> I = 1:2; W = 1:4; # set the finite domains

Model Definition

Model Initialization

The first thing we need to do is initialize our InfiniteModel and assign an appropriate optimizer that will be used to solve its transcripted variant. For our little example let's choose to use Ipopt:

julia> using InfiniteOpt, Distributions, Ipopt;

julia> model = InfiniteModel(Ipopt.Optimizer)
An InfiniteOpt Model
Feasibility problem with:
Finite Parameters: 0
Infinite Parameters: 0
Variables: 0
Derivatives: 0
Measures: 0
Optimizer model backend information:
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

Learn more about InfiniteModels and optimizers on our Infinite Models page.

Before moving on, let's go ahead make a finite parameter via @finite_parameter for $\epsilon$ since this likely a constant we'll want to update repeatedly (e.g., to determine a tradeoff curve by varying it): 

julia> @finite_parameter(model, ϵ == 10)
ϵ

Learn more about finite parameters on our Finite Parameters page.

Infinite Parameters

The next thing we need to do is identify the infinite domains our problem contains and define an infinite parameter(s) for each one via [@infinite_parameter]. For this problem we have the time domain $t \in \mathcal{D}_t$ and the random domain $\xi \in \mathcal{D}_\xi$ where $\xi \sim \mathcal{N}(\mu, \sigma^2)$:

julia> @infinite_parameter(model, t in [0, 60], num_supports = 61, 
                           derivative_method = OrthogonalCollocation(2))
t

julia> @infinite_parameter(model, ξ ~ Normal(μ, σ^2), num_supports = 10)
ξ

Notice we specify the domain/distribution the parameter depends on via in. Here we also specify the number of finite supports we desire for each parameter that will ultimately be used to reformulate and solve the problem (i.e., discretize). We also specify the derivative evaluation method associated with t that will be used evaluate the derivatives numerically. See more information about parameters on our Infinite Parameters page. Also learn more about derivative methods on our Derivative Operators page.

Variables

Now that we have an InfiniteModel and infinite parameters let's define our decision variables. First, infinite variables (ones that depend on infinite parameters) are defined via @variable with the addition of the Infinite variable type argument to specify the infinite parameters it depends on:

julia> @variable(model, x[I], Infinite(t, ξ), start = 0)
1-dimensional DenseAxisArray{GeneralVariableRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{GeneralVariableRef}:
 x[1](t, ξ)
 x[2](t, ξ)

julia> @variable(model, v[I], Infinite(t, ξ), start = 0)
1-dimensional DenseAxisArray{GeneralVariableRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{GeneralVariableRef}:
 v[1](t, ξ)
 v[2](t, ξ)

julia> @variable(model, u[I], Infinite(t), start = 0)
1-dimensional DenseAxisArray{GeneralVariableRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{GeneralVariableRef}:
 u[1](t)
 u[2](t)

julia> @variable(model, y[W] >= 0, Infinite(ξ), start = 0)
1-dimensional DenseAxisArray{GeneralVariableRef,1,...} with index sets:
    Dimension 1, 1:4
And data, a 4-element Vector{GeneralVariableRef}:
 y[1](ξ)
 y[2](ξ)
 y[3](ξ)
 y[4](ξ)

Notice that we specifying 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:

  • Finite variables: variables that do not depend on infinite parameters (defined using @variable)
  • Semi-infinite variables: infinite variables where 1 or more parameters are set a particular point (defined using @variable with the SemiInfinite variable type argument)
  • Point variables: infinite variables at a particular point (defined using @variable with the Point variable type argument).

Objective & Constraints

Now that the variables and parameters are ready to go, let's define our problem. First, we can define the objective using @objective:

julia> @objective(model, Min, integral(sum(u[i]^2 for i in I), t))
∫{t ∈ [0, 60]}[u[1](t)² + u[2](t)²]

Notice that we also employ integral to define the integral. Note that objectives must evaluate over all included infinite domains.

Now let's define the initial conditions using @constraint in combination with Restricted Variables which will restrict the domain of the variables to only be enforced at the initial time:

julia> @constraint(model, [i in I], x[i](0, ξ) == x0[i])
1-dimensional DenseAxisArray{InfOptConstraintRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{InfOptConstraintRef}:
 x[1](0, ξ) = 0, ∀ ξ ~ Normal
 x[2](0, ξ) = 0, ∀ ξ ~ Normal

julia> @constraint(model, [i in I], v[i](0, ξ) == v0[i])
1-dimensional DenseAxisArray{InfOptConstraintRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{InfOptConstraintRef}:
 v[1](0, ξ) = 0, ∀ ξ ~ Normal
 v[2](0, ξ) = 0, ∀ ξ ~ Normal

Note it is important that we include appropriate boundary conditions when using derivatives in our model. For more information please see Derivative Operators.

Next, we can add our model constraints that have derivatives using @constraint and deriv:

julia> @constraint(model, c1[i in I], deriv(x[i], t) == v[i])
1-dimensional DenseAxisArray{InfOptConstraintRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{InfOptConstraintRef}:
 c1[1] : ∂/∂t[x[1](t, ξ)] - v[1](t, ξ) = 0, ∀ t ∈ [0, 60], ξ ~ Normal
 c1[2] : ∂/∂t[x[2](t, ξ)] - v[2](t, ξ) = 0, ∀ t ∈ [0, 60], ξ ~ Normal

julia> @constraint(model, c2[i in I], ξ * deriv(v[i], t) == u[i])
1-dimensional DenseAxisArray{InfOptConstraintRef,1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{InfOptConstraintRef}:
 c2[1] : ξ*∂/∂t[v[1](t, ξ)] - u[1](t) = 0, ∀ t ∈ [0, 60], ξ ~ Normal
 c2[2] : ξ*∂/∂t[v[2](t, ξ)] - u[2](t) = 0, ∀ t ∈ [0, 60], ξ ~ Normal

Finally, we can define our last 2 constraints:

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:
    Dimension 1, 1:4
And data, a 4-element Vector{InfOptConstraintRef}:
 c3[1] : -x[1](0, ξ)² - x[2](0, ξ)² + y[1](ξ) + 2 x[1](0, ξ) + 2 x[2](0, ξ) = 2, ∀ ξ ~ Normal
 c3[2] : -x[1](25, ξ)² - x[2](25, ξ)² + y[2](ξ) + 8 x[1](25, ξ) + 6 x[2](25, ξ) = 25, ∀ ξ ~ Normal
 c3[3] : -x[1](50, ξ)² - x[2](50, ξ)² + y[3](ξ) + 12 x[1](50, ξ) = 36, ∀ ξ ~ Normal
 c3[4] : -x[1](60, ξ)² - x[2](60, ξ)² + y[4](ξ) + 2 x[1](60, ξ) + 2 x[2](60, ξ) = 2, ∀ ξ ~ Normal

julia> @constraint(model, c4, expect(sum(y[w] for w in W), ξ) <= ϵ)
c4 : 𝔼{ξ}[y[1](ξ) + y[2](ξ) + y[3](ξ) + y[4](ξ)] - ϵ ≤ 0

Notice we are able to invoke an expectation simply by calling expect.

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!:

julia> optimize!(model)

We can check the solution status via termination_status:

julia> termination_status(model)
LOCALLY_SOLVED::TerminationStatusCode = 4

Thus, our model was solved successfully! For more information please see our Optimization and Results pages.

Query the Solution

Finally, we can query a wide variety of information about our solution. Perhaps most commonly we'll want to know the objective value and the optimal primal values of decision variables. This is accomplished via objective_value and value:

julia> opt_obj = objective_value(model);

julia> u_opt = value.(u);

Note that u_opt will be multi-dimensional combination with the support values used to transcribe u(t) along the domain of t. We can query those corresponding support values via supports:

julia> u_ts = supports.(u)
1-dimensional DenseAxisArray{Vector{Tuple},1,...} with index sets:
    Dimension 1, 1:2
And data, a 2-element Vector{Vector{Tuple}}:
 [(0.0,), (1.0,), (2.0,), (3.0,), (4.0,), (5.0,), (6.0,), (7.0,), (8.0,), (9.0,)  …  (51.0,), (52.0,), (53.0,), (54.0,), (55.0,), (56.0,), (57.0,), (58.0,), (59.0,), (60.0,)]
 [(0.0,), (1.0,), (2.0,), (3.0,), (4.0,), (5.0,), (6.0,), (7.0,), (8.0,), (9.0,)  …  (51.0,), (52.0,), (53.0,), (54.0,), (55.0,), (56.0,), (57.0,), (58.0,), (59.0,), (60.0,)]

Please see the Results page for more information.

Summary Script

The example used in the sections above is summarized in the script below:

using InfiniteOpt, Distributions, Ipopt

# DEFINE THE PROBLEM CONSTANTS
μ = 1; σ = 0.2
x0 = [0, 0]; v0 = [0, 0]
p = [1 4 6 1; 1 3 0 1]; tw = [0, 25, 50, 60]
I = 1:2; W = 1:4

# INITIALIZE THE MODEL
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))
@infinite_parameter(model, ξ ~ Normal(μ, σ^2), num_supports = 10)

# INITIALIZE THE VARIABLES
@variable(model, x[I], Infinite(t, ξ), start = 0)
@variable(model, v[I], Infinite(t, ξ), start = 0)
@variable(model, u[I], Infinite(t), start = 0)
@variable(model, y[W] >= 0, Infinite(ξ), start = 0)

# SET THE OBJECTIVE
@objective(model, Min, integral(sum(u[i]^2 for i in I), t))

# SET THE INITIAL CONDITIONS
@constraint(model, [i in I], x[i](0, ξ) == x0[i])
@constraint(model, [i in I], v[i](0, ξ) == v0[i])

# SET THE PROBLEM CONSTRAINTS
@constraint(model, c1[i in I], @deriv(x[i], t) == v[i])
@constraint(model, c2[i in I], ξ * @deriv(v[i], t) == u[i])
@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), ξ) <= ϵ)

# SOLVE THE MODEL
optimize!(model)

# GET THE RESULTS
termination_status(model)
opt_obj = objective_value(model)
u_opt = value.(u)
u_ts = supports.(u)