The JuliaSim software is available for preview only. Please contact us for access, by emailing [email protected], if you are interested in evaluating JuliaSim.
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 = -k₁ * y₁ + k₃ * y₂ * y₃ du = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃ du = 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
surralg = LPCTESN(1000) sim = DEProblemSimulation(prob)
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
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).
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:
Finally, a comprehensive diagnostic report can be obtained by calling the
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: