Параллельное ансамблевое моделирование

Выполнение моделирования методом Монте-Карло, решение с предопределенным набором начальных условий, GPU-распараллеливание поиска по параметрам — все это относится к интерфейсу ансамблевого моделирования. Этот интерфейс позволяет объявить шаблон DEProblem для распараллеливания, настроить шаблон в trajectories для множества траекторий, решить каждую из них в параллельных пакетах, свести решения к конкретным ответам и вычислить суммарную статистику по результатам.

Выполнение ансамблевого моделирования

Построение задачи

Чтобы выполнить моделирование по ансамблю траекторий, определим задачу EnsembleProblem. Конструктор имеет следующий вид.

EnsembleProblem(prob::DEProblem;
                output_func = (sol, i) -> (sol, false),
                prob_func = (prob, i, repeat) -> (prob),
                reduction = (u, data, I) -> (append!(u, data), false),
                u_init = [], safetycopy = prob_func !== DEFAULT_PROB_FUNC)

output_func. Функция определяет, что будет сохранено из решения в выходной массив. По умолчанию сохраняется само решение. Выходными данными являются (out,rerun), где out — это вывод, а rerun — логическое значение, которое указывает на необходимость повторного выполнения.
prob_func. Функция, с помощью которой должна быть изменена задача. prob — задача, i — уникальный идентификатор 1:trajectories для задачи, а repeat — итерация повторения. Сначала это 1, но если rerun имело значение true, это будет 2, 3 и т. д., считая, сколько раз повторялась задача i.
reduction. Эта функция определяет способ сокращения данных в каждом пакете. По умолчанию данные (data) добавляются в u, инициализированный с помощью u_data, из пакетов. I — диапазон индексов, задающих траектории, соответствующие пакетам. Вторая часть вывода определяет, сходится ли моделирование . Если true, моделирование завершится досрочно. По умолчанию всегда имеет значение false.
u_init: Начальная форма объекта, которая обновляется на месте внутри функции reduction.
safetycopy. Определяет, вызывается ли безопасная функция deepcopy для prob перед prob_func. По умолчанию это верно для любой пользовательской функции prob_func, так как без этого изменение аргументов чего-либо в функции prob_func, например параметров или кэшей, хранящихся внутри пользовательской функции, необязательно потокобезопасно. Если вы знаете, что функция является потокобезопасной, задание этому параметру значения false может повысить производительность при использовании с потоками. Для вложенных задач, например задач SDE с пользовательскими процессами шума, функции deepcopy может оказаться недостаточно. В таких случаях используйте пользовательскую функцию prob_func.

Можно задать функцию prob_func, которая изменит задачу. Например:

function prob_func(prob, i, repeat)
    @. prob.u0 = randn() * prob.u0
    prob
end

изменяет начальное условие для всех задач на стандартное нормальное случайное число (для каждого моделирования свое случайное число). Обратите внимание, что поскольку типы задач являются неизменяемыми, она использует .=. В противном случае можно просто создать новый тип задачи:

function prob_func(prob, i, repeat)
    @. prob.u0 = u0_arr[i]
    prob
end

Если функция представляет собой ParameterizedFunction, для выполнения поиска по параметрам можно произвести аналогичные изменения в prob.f. output_func является функцией редуцирования. Ее аргументами являются сгенерированное решение и уникальный индекс для выполнения. Например, если мы хотим сохранять только 2-ю координату в конце каждого решения, можно сделать следующее.

output_func(sol, i) = (sol[end, 2], false)

Таким образом, ансамблевое моделирование будет возвращать в качестве данных массив, представляющий собой конечное значение 2-й зависимой переменной для каждого из выполнений.

Решение задачи

sim = solve(prob, alg, ensemblealg, kwargs...)

Именованные аргументы принимают аргументы для общего интерфейса решателя и передают их решателю дифференциальных уравнений. Параметр ensemblealg является необязательным и по умолчанию принимает значение EnsembleThreads(). Следует обратить внимание на следующие специальные именованные аргументы.

trajectories: количество моделирований для выполнения. Это обязательный аргумент.
batch_size: размер пакетов, на которые распространяются сокращения. Значение по умолчанию — trajectories.
pmap_batch_size: размер пакетов pmap. Значение по умолчанию — batch_size÷100 > 0 ? batch_size÷100 : 1.

Алгоритмы EnsembleAlgorithm

Выбор ансамблевого алгоритма позволяет управлять способом обработки нескольких траекторий. В настоящее время существуют следующие типы ансамблевых алгоритмов.

EnsembleSerial() — без параллелизма.
EnsembleThreads() — используется по умолчанию. Использует многопоточность. Только локальный (один компьютер, общая память) параллелизм. Максимально быстрый при быстрых траекториях.
EnsembleDistributed() — использует pmap внутренним образом. Он будет использовать столько процессоров, сколько у вас процессов Julia. Для добавления дополнительных процессов используйте addprocs(n). Дополнительные сведения см. в документации по Julia. Рекомендуется для случая, когда каждый расчет траектории выполняется «не слишком быстро» (как минимум около миллисекунды каждый?).
EnsembleSplitThreads() — использует потоковую обработку для каждого процесса, что позволяет разделить задачу на равные части nprocs(). Это нужно для решения множества быстрых траекторий на многоузловом компьютере. Рекомендуется иметь по одному процессу на каждом узле.
EnsembleGPUArray() — требует установки и using DiffEqGPU. При этом используется GPU для вычисления ансамбля с гиперпараллелизмом. Он автоматически перекомпилирует функции Julia в GPU. Производительность стандартного GPU по сравнению с 16-ядерным процессором Xeon выше в пять раз. Однако существуют ограничения, определяющие, какие функции могут выполнять автокомпиляцию подобным образом. Более подробную информацию см. в файле сведений DiffEqGPU.

Например, EnsembleThreads() вызывается следующим образом:

solve(ensembleprob, alg, EnsembleThreads(); trajectories = 1000)

Тип решения

Результирующим типом является EnsembleSimulation, который включает в себя массив решений.

Шаблон графика

Для AbstractEnsembleSimulation существует шаблон графика, в котором собраны все шаблоны для компонентных решений. При этом также передаются именованные аргументы. Полезным аргументом является linealpha, который изменяет прозрачность графиков. Дополнительным аргументом является idxs, который позволяет выбрать, какие компоненты решения строить. Например, если дифференциальное уравнение представляет собой вектор из 9 значений, то при idxs=1:2:9 будут построены решения только нечетных компонент. Еще один дополнительный аргумент — zcolors (псевдоним marker_z), который позволяет передавать zcolor для каждого ряда. Подробнее о zcolor см. в документации по рядам для Plots.jl.

Анализ ансамблевого эксперимента

Имеющиеся средства анализа позволяют формировать сводную статистику и сводные графики для EnsembleSimulation.

Чтобы воспользоваться этой функциональностью, импортируйте модуль анализа следующим образом:

using DifferentialEquations.EnsembleAnalysis

(или напрямую с помощью SciMLBase.EnsembleAnalysis).

Временные шаги и временные точки

Существует два типа построения сводных статистик. Можно формировать сводки либо по временным шагам, либо по временным точкам. Обобщение по временным шагам предполагает, что все временные шаги являются одной и той же временной точкой, т. е. интегратор использовал фиксированное значение dt или значения были сохранены с помощью saveat. Обобщение по временным точкам требует интерполяции решения.

Анализ на временном шаге или временной точке

get_timestep(sim, i) # Возвращает итератор каждого моделирования на временном шаге i
get_timepoint(sim, t) # Возвращает итератор каждого моделирования во временной точке i
componentwise_vectors_timestep(sim, i) # Возвращает вектор каждого моделирования на временном шаге i
componentwise_vectors_timepoint(sim, t) # Возвращает вектор каждого моделирования во временной точке i

Сводная статистика

Однократная статистика

Доступны следующие функции для временных шагов.

timestep_mean(sim, i) # Вычисляет среднее значение каждой компоненты на временном шаге i
timestep_median(sim, i) # Вычисляет медиану каждой компоненты на временном шаге i
timestep_quantile(sim, q, i) # Вычисляет квантиль q каждой компоненты на временном шаге i
timestep_meanvar(sim, i)  # Вычисляет среднее значение и дисперсию каждой компоненты на временном шаге i
timestep_meancov(sim, i, j) # Вычисляет среднее значение в i и j, а также ковариацию каждой компоненты
timestep_meancor(sim, i, j) # Вычисляет среднее значение в i и j, а также корреляцию каждой компоненты
timestep_weighted_meancov(sim, W, i, j) # Вычисляет среднее значение в i и j, а также взвешенную ковариацию W каждой компоненты

Доступны следующие функции для временных точек.

timepoint_mean(sim, t) # Вычисляет среднее значение каждой компоненты во временной точке t
timepoint_median(sim, t) # Вычисляет медиану каждой компоненты во временной точке t
timepoint_quantile(sim, q, t) # Вычисляет квантиль q каждой компоненты во временной точке t
timepoint_meanvar(sim, t) # Вычисляет среднее значение и дисперсию каждой компоненты во временной точке t
timepoint_meancov(sim, t1, t2) # Вычисляет среднее значение во временной точке t1 и t2, а также ковариацию каждой компоненты
timepoint_meancor(sim, t1, t2) # Вычисляет среднее значение во временной точке t1 и t2, а также корреляцию каждой компоненты
timepoint_weighted_meancov(sim, W, t1, t2) # Вычисляет среднее значение во временной точке t1 и t2, а также взвешенную ковариацию W каждой компоненты

Статистика полных временных рядов

Кроме того, для анализа полного временного ряда доступны следующие функции. Версии mean и meanvar возвращают значение DiffEqArray, которое может быть построено напрямую. meancov и meancor возвращают матрицу кортежей, где кортежи — это (mean_t1,mean_t2,cov or cor).

Доступны следующие функции для временных шагов.

timeseries_steps_mean(sim) # Вычисляет среднее значение на каждом временном шаге
timeseries_steps_median(sim) # Вычисляет медиану на каждом временном шаге
timeseries_steps_quantile(sim, q) # Вычисляет квантиль q на каждом временном шаге
timeseries_steps_meanvar(sim) # Вычисляет среднее значение и дисперсию на каждом временном шаге
timeseries_steps_meancov(sim) # Вычисляет матрицу ковариации и средние значения на каждом временном шаге
timeseries_steps_meancor(sim) # Вычисляет матрицу корреляции и средние значения на каждом временном шаге
timeseries_steps_weighted_meancov(sim) # Вычисляет матрицу взвешенной ковариации и средние значения на каждом временном шаге

Доступны следующие функции для временных точек.

timeseries_point_mean(sim, ts) # Вычисляет среднее значение в каждой временной точке в ts
timeseries_point_median(sim, ts) # Вычисляет медиану в каждой временной точке в ts
timeseries_point_quantile(sim, q, ts) # Вычисляет квантиль q в каждой временной точке в ts
timeseries_point_meanvar(sim, ts) # Вычисляет среднее значение и дисперсию в каждой временной точке в ts
timeseries_point_meancov(sim, ts) # Вычисляет матрицу ковариации и средние значения в каждой временной точке в ts
timeseries_point_meancor(sim, ts) # Вычисляет матрицу корреляции и средние значения в каждой временной точке в ts
timeseries_point_weighted_meancov(sim, ts) # Вычисляет матрицу взвешенной ковариации и средние значения в каждой временной точке в ts

EnsembleSummary

Тип EnsembleSummary включен для помощи при анализе общей сводной статистики. Имеются два конструктора:

EnsembleSummary(sim; quantiles = [0.05, 0.95])
EnsembleSummary(sim, ts; quantiles = [0.05, 0.95])

Первый создает сводку (mean,var) на каждом временном шаге. Как и в случае со сводной статистикой, здесь предполагается, что временные шаги одинаковы. Второй создает сводку (mean,var) в каждой временной точке t в ts. Для этого требуется возможность интерполяции решения. Квантиль используется для определения квантилей qlow и qhigh в каждой временной точке. По умолчанию используются квантили 5 % и 95 %.

Шаблон графика

В состав EnsembleSummary входит шаблон графика для визуализации сводной статистики. Доступны следующие дополнительные именованные аргументы.

idxs: компоненты решения, график которых будет построен. По умолчанию строятся графики всех компонент.
error_style: стиль для построения графика погрешности. Значение по умолчанию — ribbon. Другие варианты: :bars для полос погрешностей и :none, когда полосы погрешностей отсутствуют.
ci_type: По умолчанию используется значение :quantile, которое имеет квантили (qlow,qhigh), ограничения которых были определены при построении EnsembleSummary. Доверительный интервал Гаусса 1.96*(standard error of the mean) может быть задан с помощью ci_type=:SEM.

Одним из полезных аргументов является fillalpha, который управляет прозрачностью ленты вокруг среднего значения.

Пример 1. Решение уравнения ODE с различными начальными условиями

Случайные начальные условия

Проверим чувствительность линейного уравнения ODE к его начальному условию. Для этого необходимо 100 раз решить линейное уравнение ODE и построить графики траекторий. Для начала откроем несколько дополнительных процессов, чтобы таким образом распараллелить вычисления. Здесь мы будем использовать распределенный параллелизм, что означает, что необходимые функции должны быть доступны всем процессам. Это можно сделать с помощью макроса @everywhere:

using Distributed
using DifferentialEquations
using Plots

addprocs()
@everywhere using DifferentialEquations

Теперь определим линейное уравнение ODE, которое является нашей базовой задачей.

# Линейное уравнение ODE, которое начинается в точке 0.5 и решается от t=0.0 до t=1.0
prob = ODEProblem((u, p, t) -> 1.01u, 0.5, (0.0, 1.0))

Для ансамблевого моделирования нам нужно изменить начальное условие. Для этого используем функцию prob_func. Эта функция принимает базовую задачу и изменяет ее для создания новой задачи, которую в действительности решает траектория. Здесь мы возьмем базовую задачу, умножим начальное условие на rand() и используем это для вычисления траектории:

@everywhere function prob_func(prob, i, repeat)
    remake(prob, u0 = rand() * prob.u0)
end

Теперь мы строим и решаем EnsembleProblem с помощью этой базовой задачи и prob_func:

ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
sim = solve(ensemble_prob, Tsit5(), EnsembleDistributed(), trajectories = 10)

С помощью шаблона графика можно построить график того, как выглядят 10 уравнений ODE:

plot(sim, linealpha = 0.4)

Отметим, что если нужно узнать, каким было начальное условие для данной траектории, можно извлечь его из решения. sim[i] возвращает i-й объект решения. sim[i].prob — это задача, которую решала конкретная траектория, а sim[i].prob.u0 — начальное условие, используемое в i-й траектории.

Примечание. Если в задаче есть обратные вызовы, функции для condition и affect! должны быть именованными функциями (не анонимными).

Использование многопоточности

Предыдущее ансамблевое моделирование также может быть распараллелено с использованием многопоточного подхода, что позволит задействовать различные ядра одного компьютера. Поскольку память является общей для различных потоков, использовать макрос @everywhere необязательно. Та же задача может быть реализована просто следующим образом:

using DifferentialEquations
prob = ODEProblem((u, p, t) -> 1.01u, 0.5, (0.0, 1.0))
function prob_func(prob, i, repeat)
    remake(prob, u0 = rand() * prob.u0)
end
ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
sim = solve(ensemble_prob, Tsit5(), EnsembleThreads(), trajectories = 10)
using Plots;
plot(sim);

Количество используемых потоков должно быть определено вне Julia в переменной среды JULIA_NUM_THREADS (подробнее см. в документации по Julia).

Предопределенные начальные условия

Часто вам может быть уже известно, какие начальные условия необходимо использовать. Их можно задать с помощью аргумента i функции prob_func. Этот аргумент i является уникальным индексом каждой траектории. Так, если у нас есть trajectories=100, у нас есть i как некоторый индекс в 1:100, причем для каждой траектории он свой.

Поэтому если бы мы хотели использовать сетку равномерно распределенных начальных условий от 0 до 1, мы могли бы просто проиндексировать тип linspace:

initial_conditions = range(0, stop = 1, length = 100)
function prob_func(prob, i, repeat)
    remake(prob, u0 = initial_conditions[i])
end

prob_func (generic function with 1 method)

Стоит отметить, что при успешном выполнении этого кода видимого вывода не будет.

Пример 2. Решение уравнения SDE с разными параметрами

Решим то же уравнение SDE, но с изменяющимися параметрами. Создадим систему Лотки-Вольтерра с мультипликативным шумом. В системе Лотки-Вольтерра он будет дрейфовой компонентой:

function f(du, u, p, t)
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = -3 * u[2] + u[1] * u[2]
end

f (generic function with 1 method)

В качестве функции шума мы будем использовать мультипликативный шум:

function g(du, u, p, t)
    du[1] = p[3] * u[1]
    du[2] = p[4] * u[2]
end

g (generic function with 1 method)

Теперь построим уравнение SDE с помощью этих функций:

using DifferentialEquations
p = [1.5, 1.0, 0.1, 0.1]
prob = SDEProblem(f, g, [1.0, 1.0], (0.0, 10.0), p)

SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 1.0
 1.0

Это базовая задача для нашего исследования. В этом эксперименте мы хотели бы сохранить те же параметры в детерминированной компоненте, но варьировать параметры для количества шума, используя 0.3rand(2) в качестве параметров. И снова мы делаем это с помощью функции prob_func, причем здесь мы изменяем параметры в prob.p:

# `p` является глобальной переменной, ее указание было бы неустойчивым по типу.
# Использование блока let определяет небольшую локальную область, в которой мы можем
# захватить эту локальную `p`, которая нигде в этой локальной области не переопределена.
# Так она может быть устойчивой по типу.
prob_func = let p = p
    (prob, i, repeat) -> begin
        x = 0.3rand(2)
        remake(prob, p = [p[1], p[2], x[1], x[2]])
    end
end

#1 (generic function with 1 method)

Теперь решим задачу 10 раз и построим все траектории в фазовом пространстве:

ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
sim = solve(ensemble_prob, SRIW1(), trajectories = 10)
using Plots;
plot(sim, linealpha = 0.6, color = :blue, idxs = (0, 1), title = "Phase Space Plot");
plot!(sim, linealpha = 0.6, color = :red, idxs = (0, 2), title = "Phase Space Plot")

Затем мы можем обобщить эту информацию с помощью графика EnsembleSummary с границами среднего и дисперсии. Возьмем среднее значение / квантиль через каждые 0.1 единиц времени и построим график сводки:

summ = EnsembleSummary(sim, 0:0.1:10)
plot(summ, fillalpha = 0.5)

Обратите внимание, что здесь мы использовали границы квантиля, которые по умолчанию равны [0.05,0.95] в конструкторе EnsembleSummary. Можно перейти к стандартной погрешности границ среднего, используя ci_type=:SEM в шаблоне графика.

Пример 3. Использование редуцирования для остановки при достижении средством оценки пределов допуска

Здесь мы будем решать уравнение столько раз, сколько необходимо для того, чтобы стандартная погрешность среднего для последней временной точки не превышала допуска в 0.5. Поскольку нас интересует только конечная точка, можно указать функции output_func, чтобы она отбрасывала остальные данные.

function output_func(sol, i)
    last(sol), false
end

output_func (generic function with 1 method)

Функция prob_func просто рандомизирует начальное условие:

using DifferentialEquations
# Линейное уравнение ODE, которое начинается в точке 0.5 и решается от t=0.0 до t=1.0
prob = ODEProblem((u, p, t) -> 1.01u, 0.5, (0.0, 1.0))

function prob_func(prob, i, repeat)
    remake(prob, u0 = rand() * prob.u0)
end

prob_func (generic function with 1 method)

Функция редуцирования будет добавлять данные из текущего пакета в предыдущий и объявит о сходимости, если стандартная погрешность среднего будет вычислена как достаточно низкая:

using Statistics
function reduction(u, batch, I)
    u = append!(u, batch)
    finished = (var(u) / sqrt(last(I))) / mean(u) < 0.5
    u, finished
end

reduction (generic function with 1 method)

Затем можно определить и решить задачу:

prob2 = EnsembleProblem(prob, prob_func = prob_func, output_func = output_func,
                        reduction = reduction, u_init = Vector{Float64}())
sim = solve(prob2, Tsit5(), trajectories = 10000, batch_size = 20)

EnsembleSolution Solution of length 20 with uType:
Float64

Поскольку batch_size=20, это означает, что каждые 20 моделирований будет браться этот пакет, результаты будут добавляться в предыдущий пакет, будет вычисляться (var(u)/sqrt(last(I)))/mean(u) и, если значение достаточно мало, моделирование будет завершено. В этом случае моделирование завершается только после 20 раз (т. е. после вычисления первого пакета). Это может сэкономить массу времени.

Помимо экономии времени за счет проверки сходимости, можно экономить память за счет сокращения между пакетами. Например, предположим, что нас снова интересует только среднее значение в конце. Вместо сохранения решения в конце для каждой траектории можно сохранить текущее суммирование конечных точек:

function reduction(u, batch, I)
    u + sum(batch), false
end
prob2 = EnsembleProblem(prob, prob_func = prob_func, output_func = output_func,
                        reduction = reduction, u_init = 0.0)
sim2 = solve(prob2, Tsit5(), trajectories = 100, batch_size = 20)

EnsembleSolution Solution of length 1 with uType:
Float64

конечные точки будут суммироваться после каждых 20 решений с сохранением текущей суммы. В конечном результате sim2.u будет просто числом, а sim2.u/100, соответственно, средним значением.

Пример 4. Использование средств анализа

В этом примере мы покажем, как проанализировать EnsembleSolution. Сначала сгенерируем эксперимент, проводимый методом Монте-Карло, с 10 решениями. Для решения задачи будем использовать систему 4x2 линейных стохастических дифференциальных уравнений:

function f(du, u, p, t)
    for i in 1:length(u)
        du[i] = 1.01 * u[i]
    end
end
function σ(du, u, p, t)
    for i in 1:length(u)
        du[i] = 0.87 * u[i]
    end
end
using DifferentialEquations
prob = SDEProblem(f, σ, ones(4, 2) / 2, (0.0, 1.0)) #prob_sde_2Dlinear

SDEProblem with uType Matrix{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 4×2 Matrix{Float64}:
 0.5  0.5
 0.5  0.5
 0.5  0.5
 0.5  0.5

Чтобы решить задачу 10 раз, используем конструктор EnsembleProblem и решаем с помощью trajectories=10. Поскольку мы хотим сравнивать значения на временных шагах, необходимо убедиться, что все шаги выполняются в одно и то же время. Таким образом, мы задаем adaptive=false и явным образом задаем dt.

prob2 = EnsembleProblem(prob)
sim = solve(prob2, SRIW1(), dt = 1 // 2^(3), trajectories = 10, adaptive = false)

EnsembleSolution Solution of length 10 with uType:
RODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 3, Float64, Matrix{Float64}, Matrix{Float64}, Vector{Matrix{Float64}}, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Matrix{Float64}, Matrix{Float64}}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, Nothing, SDEFunction{true, SciMLBase.FullSpecialize, typeof(Main.f), typeof(Main.σ), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, typeof(Main.σ), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Nothing}, SRIW1, StochasticDiffEq.LinearInterpolationData{Vector{Matrix{Float64}}, Vector{Float64}}, DiffEqBase.Stats, Nothing}

Обратите внимание, что если не выполнять вычисления timeseries_steps, этот код совместим с адаптивной дискретизацией по времени. Использование адаптивности обычно более эффективно.

Вычислить среднее значение и дисперсию на 3-м временном шаге можно следующим образом:

using DifferentialEquations.EnsembleAnalysis
m, v = timestep_meanvar(sim, 3)

([0.6099069274459361 0.6293846479505159; 0.7803438431424009 0.67560338640098; 0.6974586199618295 0.7006333021379738; 0.5779647058303172 0.5255352135254947], [0.2149409863166183 0.05237026372964025; 0.1202443444995385 0.06979241803596177; 0.08530534367199581 0.05422988023029949; 0.08075663715397942 0.13447477859818402])

или вычислить среднее значение и дисперсию в момент t=0.5 можно следующим образом:

m, v = timepoint_meanvar(sim, 0.5)

([0.9178257488493667 0.8132463718308369; 1.0743270195430528 0.7161778666209464; 0.9992894513675031 0.7635438171389558; 0.7837552618741324 0.5674800632634518], [1.0038177691931804 0.2948757550863388; 0.9903081093409329 0.09586256947321561; 0.25749609542416096 0.16860903591349136; 0.4086260556118997 0.16176458400406693])

Мы можем получить ряд для среднего и дисперсии на каждом временном шаге следующим образом:

m_series, v_series = timeseries_steps_meanvar(sim)

(DiffEqArray{Float64, 3, Vector{Matrix{Float64}}, Vector{Float64}, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing, Nothing}([[0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.4590354611659848 0.504878430663902; 0.6954195160668808 0.5617844570000997; 0.5867663523717925 0.5974188501458414; 0.5804463932463718 0.5046159839117782], [0.6099069274459361 0.6293846479505159; 0.7803438431424009 0.67560338640098; 0.6974586199618295 0.7006333021379738; 0.5779647058303172 0.5255352135254947], [0.9187727430579399 0.689009442719635; 0.7832524703057391 0.7326082518826428; 0.9093270383584899 0.72489447791006; 0.5723191172675123 0.6068245366328354], [0.9178257488493667 0.8132463718308369; 1.0743270195430528 0.7161778666209464; 0.9992894513675031 0.7635438171389558; 0.7837552618741324 0.5674800632634518], [1.2467442660230896 0.863850863077547; 1.1988749390095488 0.7678314756764205; 1.0881971920274691 0.9296194326126226; 0.9135928604552026 0.7026338373279828], [1.1736131974330424 0.9170922685242343; 1.1727550147907322 0.8374599410010588; 1.6630933123358447 1.090938710439789; 0.8190263048058021 0.7646137728229288], [1.070496352588129 1.0255995011992645; 1.2932363946706358 0.7957923897944246; 1.904929063164664 1.1841910279521317; 1.0341028626551678 0.9246392186103818], [1.2335892370833763 1.1053371297132013; 1.3314747743008253 1.0855798098573368; 2.436870287174397 1.2919855060696634; 1.003301357900119 1.0911512345184846]], [0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0], SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}(nothing, nothing, nothing), nothing, nothing), DiffEqArray{Float64, 3, Vector{Matrix{Float64}}, Vector{Float64}, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Nothing, Nothing}([[0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0], [0.019862557126485894 0.021524781838826935; 0.022033767766715676 0.020190317050346716; 0.03252217467125499 0.015847864680242987; 0.06022566566016834 0.0762586297609605], [0.2149409863166183 0.05237026372964025; 0.1202443444995385 0.06979241803596177; 0.08530534367199581 0.05422988023029949; 0.08075663715397942 0.13447477859818402], [0.8647026153099191 0.20005349320888435; 0.25102711462222094 0.07281268774983589; 0.03446558633584287 0.1016181399487518; 0.10516258137383064 0.17862832026672953], [1.0038177691931804 0.2948757550863388; 0.9903081093409329 0.09586256947321561; 0.25749609542416096 0.16860903591349136; 0.4086260556118997 0.16176458400406693], [2.2736410323024416 0.4386617217264339; 1.130765669303827 0.11516378411484483; 0.34353554631984734 0.17070043320573197; 0.47976950530064916 0.29567681479999597], [1.384636147938224 0.40046933841383014; 0.5521750288816459 0.23874620874394076; 2.271825216437087 0.25625614064929403; 0.4451491596230611 0.5167281877127613], [0.5017591670928092 0.539640360783523; 0.7338733483992765 0.23597454211335694; 2.5841101167153964 0.27191171560246763; 1.4359785802694278 1.0550981578063539], [0.9080154463331953 0.4671970117761702; 0.8137468772900661 0.9886787336033822; 5.796535069183165 0.25411287974136376; 1.503452572277191 1.3437110895162936]], [0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0], SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}(nothing, nothing, nothing), nothing, nothing))

или при выбранных значениях t:

ts = 0:0.1:1
m_series = timeseries_point_mean(sim, ts)

t: 0.0:0.1:1.0
u: 11-element Vector{Matrix{Float64}}:
 [0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5]
 [0.46722836893278785 0.5039027445311217; 0.6563356128535046 0.5494275656000798; 0.5694130818974341 0.577935080116673; 0.5643571145970976 0.5036927871294226]
 [0.5495583409339557 0.5795821610358703; 0.7463741123121928 0.6300758146406279; 0.6531817129258146 0.6593475213411207; 0.5789573807967392 0.517167521680008]
 [0.7334532536907375 0.6532345658581635; 0.7815072940077362 0.6984053325936451; 0.7822059873204936 0.7103377724468084; 0.5757064704051953 0.558050942768431]
 [0.9185833442162255 0.7138568285418754; 0.8414673801532018 0.7293221748303035; 0.9273195209602925 0.7326243457558392; 0.6146063461888363 0.5989556419589588]
 [0.9178257488493667 0.8132463718308369; 1.0743270195430528 0.7161778666209464; 0.9992894513675032 0.7635438171389557; 0.7837552618741322 0.5674800632634518]
 [1.180960562588345 0.853729964828205; 1.1739653551162497 0.7575007538653258; 1.0704156438954757 0.8964043095178893; 0.8876253407389886 0.6756030825150766]
 [1.2028656248690615 0.8957957063455593; 1.183202984478259 0.8096085548712034; 1.4331348642124941 1.0264109993089225; 0.8568529270655623 0.7398217986249505]
 [1.1323664594950773 0.9604951615942465; 1.2209475667426937 0.8207929205184051; 1.7598276126673724 1.1282396374447259; 0.9050569279455486 0.8286239511379101]
 [1.1031149294871785 1.041547026902052; 1.3008840705966738 0.8537498738070071; 2.0113173079666105 1.205749923575638; 1.0279425617041582 0.9579416217920025]
 [1.2335892370833763 1.105337129713201; 1.331474774300825 1.0855798098573368; 2.4368702871743966 1.2919855060696634; 1.0033013579001193 1.0911512345184848]

Заметим, что эти ряды средних и дисперсий можно построить напрямую. Аналогичным образом можно вычислить ковариационные матрицы:

timeseries_steps_meancov(sim) # Использование временных шагов с предположением фиксированного значения dt
timeseries_point_meancov(sim, 0:(1 // 2^(3)):1, 0:(1 // 2^(3)):1) # Использование временных точек, интерполяция

9×9 Matrix{Tuple{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}}:
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])                                          …  ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.459035 0.504878; 0.69542 0.561784; 0.586766 0.597419; 0.580446 0.504616], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])      ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [0.459035 0.504878; 0.69542 0.561784; 0.586766 0.597419; 0.580446 0.504616], [0.118775 0.0154271; 0.114333 0.0434085; -0.100862 0.0026158; 0.212745 0.228353])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.609907 0.629385; 0.780344 0.675603; 0.697459 0.700633; 0.577965 0.525535], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])     ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [0.609907 0.629385; 0.780344 0.675603; 0.697459 0.700633; 0.577965 0.525535], [0.406312 0.0707691; 0.24954 0.121217; -0.266544 0.0636461; 0.270075 0.386317])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.918773 0.689009; 0.783252 0.732608; 0.909327 0.724894; 0.572319 0.606825], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])     ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [0.918773 0.689009; 0.783252 0.732608; 0.909327 0.724894; 0.572319 0.606825], [0.797289 0.147197; 0.426011 0.214776; 0.155782 0.0477489; 0.231206 0.442482])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [0.917826 0.813246; 1.07433 0.716178; 0.999289 0.763544; 0.783755 0.56748], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])       ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [0.917826 0.813246; 1.07433 0.716178; 0.999289 0.763544; 0.783755 0.56748], [0.849325 0.219076; 0.797854 0.227508; 0.857093 0.10425; 0.701703 0.406921])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [1.24674 0.863851; 1.19887 0.767831; 1.0882 0.929619; 0.913593 0.702634], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])      …  ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [1.24674 0.863851; 1.19887 0.767831; 1.0882 0.929619; 0.913593 0.702634], [1.23275 0.3333; 0.867637 0.255481; 1.18925 0.124019; 0.76114 0.597915])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [1.17361 0.917092; 1.17276 0.83746; 1.66309 1.09094; 0.819026 0.764614], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])          ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [1.17361 0.917092; 1.17276 0.83746; 1.66309 1.09094; 0.819026 0.764614], [1.01961 0.408894; 0.593937 0.393172; 3.41236 0.125985; 0.743009 0.624961])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [1.0705 1.0256; 1.29324 0.795792; 1.90493 1.18419; 1.0341 0.924639], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])              ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [1.0705 1.0256; 1.29324 0.795792; 1.90493 1.18419; 1.0341 0.924639], [0.517254 0.458208; 0.767004 0.4735; 3.40336 0.182188; 1.43925 0.981698])
 ([0.5 0.5; 0.5 0.5; 0.5 0.5; 0.5 0.5], [1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0])              ([1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [1.23359 1.10534; 1.33147 1.08558; 2.43687 1.29199; 1.0033 1.09115], [0.908015 0.467197; 0.813747 0.988679; 5.79654 0.254113; 1.50345 1.34371])

Для общего анализа можно создать тип EnsembleSummary.

summ = EnsembleSummary(sim)

EnsembleSolution Solution of length 9 with uType:
Float64

будет суммироваться на каждом временном шаге, тогда как

summ = EnsembleSummary(sim, 0.0:0.1:1.0)

EnsembleSolution Solution of length 11 with uType:
Float64

будет суммироваться во временных точках 0.1 с использованием интерполяций. Для визуализации результатов можно построить график. Поскольку в дифференциальном уравнении 8 компонент, это может привести к путанице, поэтому построим график только для 3-й компоненты:

using Plots;
plot(summ; idxs = 3);

Можно изменить ленты на полосы погрешностей и построить график двух разных индексов:

plot(summ; idxs = (3, 5), error_style = :bars)

Или мы можем просто построить график среднего значения каждой компоненты с течением времени:

plot(summ; error_style = :none)