Планирование экспериментов

Это руководство было изначально предоставлено Арпитом Бхатией (Arpit Bhatia) и Крисом Кои (Chris Coey).

В этом руководстве рассматриваются примеры планирования экспериментов (D-оптимального, A-оптимального и E-оптимального) из раздела 7.5 работы Boyd2004.

В руководстве используются следующие пакеты.

using JuMP
import SCS
import LinearAlgebra
import Random

В этом руководстве используются множества из MathOptInterface. JuMP по умолчанию экспортирует символ MOI как псевдоним для пакета MathOptInterface.jl. Мы рекомендуем сделать это более явным в вашем коде, добавив следующие строки:

import MathOptInterface as MOI

Зададим начальное значение, чтобы случайные числа были повторяемыми:

Random.seed!(1234)

Random.TaskLocalRNG()

Ослабленная задача планирования экспериментов

В базовом виде задача планирования экспериментов заключается в следующем.

Для данного набора возможных вариантов экспериментов и общего количества экспериментов , которые необходимо провести, нужно выбрать количество экспериментов каждого типа, то есть , чтобы обеспечить небольшую (в некотором смысле) ковариацию погрешности .

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

Базовая задача планирования экспериментов может оказаться сложной комбинаторной задачей, когда , то есть общее количество экспериментов, сопоставимо с , так как в этом случае все — это небольшие целые числа.

Однако в случае, когда велико по сравнению с , хорошее приближенное решение можно найти, проигнорировав или ослабив ограничение, согласно которому — это целые числа.

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

Вектор удовлетворяет условию . Кроме того, каждое является целым числом, кратным . Если проигнорировать это последнее ограничение, мы придем к следующей задаче:

Для задачи планирования экспериментов, которая представляет собой задачу векторной оптимизации на положительно полуопределенном конусе, было предложено несколько скаляризаций.

q = 4 # размерность оценочного пространства
p = 8 # количество экспериментальных векторов
n_max = 3 # верхняя граница лямбды
n = 12

V = randn(q, p)

eye = Matrix{Float64}(LinearAlgebra.I, q, q);

A-оптимальное планирование

При A-оптимальном планировании экспериментов мы минимизируем , след ковариационной матрицы. Эта целевая функция представляет собой просто среднее нормы квадрата погрешности:

Задача A-оптимального планирования экспериментов в форме SDP имеет следующий вид:

aOpt = Model(SCS.Optimizer)
set_silent(aOpt)
@variable(aOpt, np[1:p], lower_bound = 0, upper_bound = n_max)
@variable(aOpt, u[1:q], lower_bound = 0)
@constraint(aOpt, sum(np) <= n)
for i in 1:q
    matrix = [
        V*LinearAlgebra.diagm(0 => np ./ n)*V' eye[:, i]
        eye[i, :]' u[i]
    ]
    @constraint(aOpt, matrix >= 0, PSDCone())
end
@objective(aOpt, Min, sum(u))
optimize!(aOpt)
@assert is_solved_and_feasible(aOpt)
objective_value(aOpt)

5.103223362935165

value.(np)

8-element Vector{Float64}:
 2.949521135024129
 1.7790921964770634
 8.104871349218623e-6
 4.706163791527242e-6
 2.10826698418427
 1.698346379809937
 1.2635010842547454
 2.2012292330175063

E-оптимальное планирование

При -оптимальном планировании минимизируется норма ковариационной матрицы погрешностей, то есть максимальное собственное значение .

Так как диаметр (двойная длина самой длинной полуоси) доверительного эллипсоида пропорционален , минимизацию можно геометрически интерпретировать как минимизацию диаметра доверительного эллипсоида.

E-оптимальное планирование можно также интерпретировать как минимизацию максимальной дисперсии по всем с . Задача E-оптимального планирования экспериментов в форме SDP имеет следующий вид:

eOpt = Model(SCS.Optimizer)
set_silent(eOpt)
@variable(eOpt, 0 <= np[1:p] <= n_max)
@variable(eOpt, t)
@constraint(
    eOpt,
    V * LinearAlgebra.diagm(0 => np ./ n) * V' - (t .* eye) >= 0,
    PSDCone(),
)
@constraint(eOpt, sum(np) <= n)
@objective(eOpt, Max, t)
optimize!(eOpt)
@assert is_solved_and_feasible(eOpt)
objective_value(eOpt)

0.4353863856695074

value.(np)

8-element Vector{Float64}:
  2.999998678509857
  2.7528101637391016
 -3.609558428104128e-7
 -1.487238376393808e-6
  2.1818462791653928
  2.3252635463929763
  0.21896764670154034
  1.521116305115523

D-оптимальное планирование

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

dOpt = Model(SCS.Optimizer)
set_silent(dOpt)
@variable(dOpt, np[1:p], lower_bound = 0, upper_bound = n_max)
@variable(dOpt, t)
@objective(dOpt, Max, t)
@constraint(dOpt, sum(np) <= n)
E = V * LinearAlgebra.diagm(0 => np ./ n) * V'
@constraint(dOpt, [t; 1; triangle_vec(E)] in MOI.LogDetConeTriangle(q))
optimize!(dOpt)
@assert is_solved_and_feasible(dOpt)
objective_value(dOpt)

0.30845393467249654

value.(np)

8-element Vector{Float64}:
 0.4275824266205948
 2.9100163168941116
 2.916872474296307e-6
 2.6291979913959343e-6
 2.9157806299833644
 2.673375664593826
 2.7353950122202764
 0.3378388086254706