Generating surrogates

The JuliaSim software is available for preview only. Please contact us for access, by emailing [email protected], if you are interested in evaluating JuliaSim.

Introduction

In order to demonstrate how JuliaSim.jl could be used in conjunction with OrdinaryDiffEq.jl, we construct a concise example, namely, the so-called ROBER problem, which consists of a stiff system of 3 non-linear ordinary differential equations (ODEs).

First, we define the ODE in question:

using QuasiMonteCarlo, OrdinaryDiffEq, ModelingToolkit, Surrogates, JuliaSim
function rober(du, u, p, t)
    y₁, y₂, y₃ = u # initial vector
    k₁, k₂, k₃ = p # rate constants
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] =  k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] =  k₂ * y₂^2
    nothing
end

tstop = 1e4
p = [0.04, 3e7, 1e4]
u0 = [1.0, 0.0, 0.0] # initial condition
tspan = (0.0, tstop)

prob = ODEProblem(rober, u0, tspan, p)
sol = solve(prob, Rosenbrock23())

This creates an ODEProblem, which can be simulated for a set of parameters. Now we need to specify a space of parameters over which we want to train our surrogate:

param_space = [(0.036, 0.044), (2.7e7, 3.3e7), (0.9e4, 1.1e4)]

Next, we specify the surrogate algorithm and cast the ODEProblem to DEProblemSimulation:

surralg = LPCTESN(1000)
sim = DEProblemSimulation(prob)

The abbreviation LPCTESN stands for 'Linear Projection Continuous-Time Echo State Network' (consult this paper for more information). The integer 1000 specifies reservoir size. It should always be provided by the user.

Finally, we call surrogatize and simulate_ensemble as follows:

odesurrogate = JuliaSim.surrogatize(
    sim,
    param_space,
    surralg,
    100; # 'n_sample_pts'
    ensemble_kwargs = (;),
    component_name = :robertson_surrogate,
    verbose=true)
odesimresults = JuliaSim.simulate_ensemble(sim, param_space, 100)

Note that the number of sets to be sampled from the parameter space (n_sample_pts) has to be provided by the user.

As can be inferred from this example, JuliaSim is a very intuitive tool allowing for high composability within the SciML ecosystem (and beyond). It ought to be emphasized that this simple example does not fully reflect JuliaSim's capabilities (as the paragraphs above clearly demonstrate).

Plotting functions

Numerical results would be incomplete without proper visualization and summaries. JuliaSim provides the necessary plotting functionalities for an easier analysis of the results. In its simplest form, the user may just call:

plot(odesimresults; ns = 1:2, output_dim = 2) # ns stands for 'samples to be plotted'

to obtain a visualization of the ROBER problem (see above).

The user may also wish to inspect the training statistics of the CTESN:

plot(odetrainstats, odesurrogate; ns = 1:2, output_dim = 2, log_scale = true)

Below is a sample plot: training plot

Finally, a comprehensive diagnostic report can be obtained by calling the weave_diagnostics function:

JuliaSim.weave_diagnostics(
    surrogate,
    "DCPM Temperature Diagnostic Report",  # Heading of the Report
    "output_fmu.pdf";                      # Path to save the Generated Report
    doctype = "md2pdf",                    # Any output format from markdown works.
    log_scale = false,
    include_errors = [:pointwise_error, :curve_distance],
    use_absolute_err = false,
    visualize_n_samples = 1
)

Here is a part of a sample diagnostic report of the ROBER problem:

rober diagnostic report

rober plot