Общие сведения. При построении по экспериментальным данным математических моделей исследуемых объектов, описываемых уравнением вида

При построении по экспериментальным данным математических моделей исследуемых объектов, описываемых уравнением вида

y=b^Tf(x) (5.1)

где f(х) — заданные функции: f(x)={f₀(x), f₁(x),…,f_d(x)}^T; d — число параметров уравнения (5. 1); х = { х ₁, х ₂,…, х_k } — вектор независимых управляемых переменных (факторов), обычно требуется составить программу проведения эксперимента, удовлетворяющую условиям выполнения некоторого критерия оптимальности. В планировании экспе­римента широко используется критерий D-оптимальности, позволяющий охватить широкий круг экспериментальных задач: нелинейных, последовательного планирования, с произвольной областью варьирования независимых переменных, с неодинаковой точностью опытов при различных условиях проведения эксперимента.

Пусть задан точный план эксперимента

x={x⁽¹⁾, x⁽²⁾,…,x^(N)}^T

где x⁽ⁱ⁾- i-я точка факторного пространства: x⁽ⁱ⁾=(x₁⁽ⁱ⁾, x₂⁽ⁱ⁾,…, x_k⁽ⁱ⁾); N — число точек плана.

Матрица значений функций независимых переменных в точках плана имеет вид

.

Так как часть точек плана может повториться, то его можно представить в виде

,

где x^(j) — точки, в которых сосредоточен план X, j = 1,r (спектр плана); h_j — число повторных наблюдений в j -й точке плана, причем

Обозначим l_j=h_j/N — частота j -й точки плана. Тогда соответствующий ему план, заданный в виде

называется нормированным планом. Если частоты l_j могут принимать любые значения в интервале [0,1] при условии ål_j== 1, то план L называется непрерывным.

Например, если для двух переменных х₁, х₂ точный план

-1 -1

+1 -1

-1 +1

Х= -1 +1

0 0

0 0

то с учетом кратности последних точек соответствующий ему непрерывный план будет

x⁽¹⁾, x⁽²⁾, x⁽³⁾, x⁽⁴⁾, x⁽⁵⁾

L=

1/6, 1/6, 1/6, 1/6, 2/6,

Критерий D-оптимальности требует такого выбора плана X *, содержащего N опытов, при котором определитель дисперсионной матрицы С (X) минимален т.е.

½ С (X*)½=min½ С (X)½=min½ (F^TF)^-1 ½ xÎW_x

или, соответственно, максимален определитель информационной матри­цы М(Х):

½M (X*)½=max½M (X)½=max½ F^TF ½ xÎW_x

где X* — оптимальный план в смысле критерия D-оптимальности; W _х — область изменения параметров плана X.

План L * будет непрерывным D-оптимальным планом, если он минимизирует на множестве всех непрерывных планов в области W_x величину определителя дисперсионной матрицы или, соответственно, максимизирует определитель информационной матрицы М(Х).

Из теории планирования эксперимента известно, что в случае равноточных наблюдений процедура построения непрерывного D-оптимального плана сводится к выполнению рекуррентных операций, определяемых уравнениями

M(t+1)=M(t)+f(x*)f^T(x*) (5.2)

f^T(x*)C(t)f(x*)=max f^T(x*)C(t)f(x*) (5.3)

где C(t)M(t) — дисперсионная и информационная матрицы соот­ветственно на t -м шаге.

Для упрощения вычислений целесообразно разделить процедуры получения координат точек (спектра) D-оптимального плана и определения частот повторения. Выявление точек, в которых концентрируется D-оптимальный план, можно осуществить за сравнительно небольшое количество циклов по рекуррентным соотношениям (5.2), (5,3), тогда как точное определение частот в каждой точке требует большего числа циклов. Такое разбиение процедуры получения D-оптимальных планов обусловлено тем, что при определении частоты повторения наблюдений нет необходимости искать глобальный максимума f(x)^t C(t)f(x} по всему пространству, так как заранее известно, что он будет иметь место в одной из точек, найденных на первом этапе.

Для получения спектра непрерывного D-оптимального плана следует провести следующие операции:

а) выбрать произвольный начальный невырожденный план для числа наблюдений Nо (r <=Nо:

x⁽¹⁾, x⁽²⁾,…, x^(r)

L⁽⁰⁾=

l₁, l₂, …, l_r

и определить информационную матрицу

M(0)=ål_if(x⁽ⁱ⁾)f^T(xⁱ)

б) по уравнению (5.3) определить точку х *, в которой квадратичная форма f^т(x)С(0)f(х) имеет глобальный максимум на множестве W х. Поиск глобального максимума может быть основан на многократном применении локального поиска из разных точек пространства W х и последующем выборе максимального из всех значений локальных максимумов. Поиск локальных максимумов следует начинать с точек начального плана. Для обнаружения возможных максимумов, не предусмотренных начальным планом, локальный поиск следует осуществлять также из ряда случайных точек, координаты которых можно получить с помощью генератора случайных чисел, равномерно распределенных на области W х;

в) после определения точки глобального максимума х * скорректировать матрицу L⁽⁰⁾ по уравнению (5.2), т.е. получить план

x⁽¹⁾ , x⁽²⁾, …, x^(r), x*

L⁽¹⁾=

(l-a₀)l₁, (l-a₀)l₂, …, (l-a₀)l_r, a₀

Операции «б», «в» повторить с заменой плана L⁽⁰⁾ на L⁽¹⁾ и т.д. до выполнения останова в соответствии с выбранным правилом.

Практика показывает, что обычно можно производить останов по процедуре (5.2), (5.3), когда число циклов будет в два-три раза больше максимального числа точек, в которых концентрируется D-оптимальный план. В общем случае количество вычислений зависит от выбранных начальных приближений и будет тем меньше, чем ближе М(О) к информационной матрице D-оптимального плана.

Для определения частоты повторения наблюдений в каждой точке необходимо:

а) выбрать начальный план L⁽⁰⁾, включающий по одному разу все точки, которые были определены на первом этапе:

x⁽¹⁾, x⁽²⁾,…, x^(r)

L⁽⁰⁾=

1/No, 1/No, 1/No

где No —число наблюдений (в начальном плане N ₀ = r ); r — число точек спектра D-оптимального плана;

б) на основании соотношения (5.3) определить точку x⁽ⁱ⁾* спектра плана, в которой квадратичная форма f^т(х^(i)*)C(0)f(x^(i)*) (i=1, r) больше, чем в остальных точках спектра. В том случае, когда получаются одинаковые наибольшие значения в нескольких точках спектра оптимального плана, выбирается любая из них;

в) скорректировать информационную матрицу М по уравнению (5.2), в результате чего получается план

где а₀ = 1 / N и k — номер точки спектра, в которой добавляется еще одно наблюдение к плану L ⁽⁰⁾;

— операции «б»и«в» повторить с заменой плана L⁽⁰⁾ на L⁽¹⁾ и числа наблюдений N ₀ на N ₁.

Другой подход к определению частот связан с подсчетом числа попаданий глобального максимума квадратичной формы f^т(х^(i)*)C(0)f(x^(i)*) (i=1, r) в каждую точку спектра х⁽ⁱ⁾ D- оптимального плана:

l_i = (m_i+l)/(N₀ +S),

где т i — число попаданий в точку х^(i)* спектра плана; S — число циклов на этапе определения частот.

Из процедуры построения непрерывных D-оптимальных планов видно, что вычислительные затраты на каждом шаге планирования складываются главным образом из времени, необходимого для поиска гло­бального максимума квадратичной формы, и времени обращения скорректированной информационной матрицы.

Дисперсионная матрица С (t + 1) на (t + 1)-м шаге рекуррентной процедуры может быть получена на основании известной после t-го этапа информационной матрицы по формуле

C(t+1)=[M(t)+f(x)f^T(x)]^-1

где х — точка, добавляемая в план L^(t).

Для упрощения вычислений можно воспользоваться известной формулой обращения матриц

(5.4)

где А — квадратичная матрица размерности п х п; U — вектор-столбец размерности п. С учетом того что добавление точки в исходный план происходит c некоторым весом а, выражение (5.4) приводит к следующей формуле для определения элементов обратной матрицы С (t + 1) по известной матрице С(t):

Обозначим через Q (х, L) функцию Q(x, L)=f^T(x)C(L^(t))f(x). Из теории планирования эксперимента известно, что непрерывный план L тогда и только тогда оптимален, когда

max Q(x, L) = k+1,

x Î W_x

где (k+ 1) — число оцениваемых параметров.

Это положение можно использовать для оценки относительного отличия получаемого плана L от D-оптимального с помощью формулы

С другой стороны, этой формулой можно воспользоваться для останова процедуры построения D-оптимального плана при достижении некоторого наперед заданного, достаточно малого положительного значения d.

Полученный непрерывный план может быть использован для постро­ения точного D-оптимального плана при заданном числе опытов N. Решение данной задачи зависит от соотношений числа опытов в точном плане N и числа точек непрерывного плана r, а также от соотношений максимальной l_mах и минимальной l_min частот точек непрерывного плана. При этом возможны следующие ситуации:

а) N = т r, l _mах = l _min и т — целое число. Так как l _mах = l _min, то

l₁=l₂=…=l_r=1/r

Точный план, определяемый при этих условиях с помощью непрерыв­ного D-оптимального плана, является D-оптимальным планом. Количество наблюдений в точке х⁽ⁱ⁾ этого плана

h_i=nl_i=N1/r=mr1/r=m

б) r <<N. В этом случае можно ожидать, что с помощью непрерыв­ного плана получится точный план, достаточно близкий к D-оптимальному. Число наблюдений hi в точке х⁽ⁱ⁾ определяется округлением произведения N1i до ближайшего целого числа;

в) k + 1 <= N <r. В этом случае трудно получить однозначное решение.

Варианты задания приведены в табл.5.1.

Таблица 5.1

Номер бригады	Вид математической модели объекта
	y=bo+b1x1+b2x12
	y=bo+b1x1+b2x2
	y=bo+b1x1+ b2x2+b12x1x2
	y=bo+b1x1+ b2x2+b11x12
	y=bo+b1x1+ b2x2+ b12x1x2+b11x12
	y=bo+b1x1+ b2x2+ b12x1x2+b11x12+b22x22
	y=bo+b11x12+b22x22

Построить D-оптимальные планы для расположения точек, представленного на рис.5.1.

Содержание отчета

Отчет должен содержать:

—задание;

— блок-схему процедуры вычисления непрерывного D-оптимального

плана; ооауг

— программу для вычисления D-оптимального плана на JBM;

— полученный непрерывный и точный D-оптимальные планы;

— выводы по работе.

Контрольные вопросы

1. Какой план называется D-оптимальным?

2 Чем отличается непрерывный план от точного?

3. Чем вызвано разбиение процедуры построения D-оптимального

плана на два этапа?

4. Как определяется глобальный максимум квадратичной формы f^T(x)Cf(x)?

5. Исходя из каких условий выполняется останов вычислительных процедур каждого этапа построения D-оптимального плана?

6. Каким образом проверяется близость полученного плана к D-оптимальному?

7. Как получить на основе непрерывного точный D-оптимальный план?

Приложение 1

ЗНАЧЕНИЯ F-РАСПРЕДЕЛЕНИЯ ПРИ УРОВНЕ ЗНАЧИМОСТИ а = 0,05