Планирование активного эксперимента при поиске оптимальных условий

Теоретическая часть

Одной из центральных задач планирования эксперимента является задача поиска экстремума функции отклика [1] — например, при оптимизации производственных и других процессов.

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

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

При решении задач планирования эксперимента наибольшее распространение получили алгоритмы поиска, использующие градиентные методы. Их особенность состоит в том, что движение при поиске происходит в направлении не самого градиента, который неизвестен, а его оценки. Если функция отклика априорно выражается зависимостью $$y = f(X_{1}, X_{2}, … , X_{k})$$, то оценка градиента $$grad\mbox{ }f(X_{1}, X_{2}, … , X_{k})$$ в точке факторного пространства при этом находится по результатам измерений, проводимым в ее окрестности. Из класса градиентных методов практическое применение получил метод, разработанный Боксом и Уильсоном. В нем используется метод крутого восхождения в сочетании с последовательно планируемым факторным экспериментом для нахождения оценки градиента [1]. В общем виде метод Бокса и Уильсона состоит в следующем:

построение факторного эксперимента в окрестности некоторой точки;

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

крутое восхождение в направлении оценки градиента;

нахождение оценки максимума (минимума) функции отклика по этому направлению;

анализ оценки максимума функции отклика и принятие решения о дальнейшем восхождении.

1. Оценивание градиента

Рассмотрим оценивание градиента функции отклика, когда факторы имеют два уровня. Пусть функция отклика

$$y = f(X_1, X_2, … , X_k)$$

определена в области $$G\subset R^k$$.

Используем произвольную точку факторного пространства $$X^0=(x_1^0,X_2^0,...,X_k^0)$$ как центр плана. Верхний $$X_{i2}^0$$ и нижний $$X_{i1}^0$$ уровни выбираем симметричными относительно центра плана, чтобы перейти к кодированным переменным (факторам):

$$x_i=\frac{X_i-X_i^0}{S_i^0},\mbox{ }i=1,2,...,k,$$

где $$S_i^0=(X_{i2}^0-X_{i1}^0)/2$$ — интервал варьирования.

Тогда функцию отклика (10.1) можно выразить через кодированные переменные:

$$y = f_1(х_1, х_2, … , х_k).$$

Под задачей оценивания градиента теперь понимается определение оценки градиента функции отклика (10.2) в точке $$x^0=(x_1^0,x_2^0,...,x_k^0)$$, где $$x_i^0$$ $$(i=1,2,...,k)$$, т. е. в нулевой точке кодированных переменных. Принимаем, что в окрестности точки $$(x_1^0,x_2^0,...,x_k^0)$$ функция отклика (10.2) допускает разложение по формуле Тейлора вида

$$y=f_1(x_1,x_2,...,x_k)=f_1(x_1^0,x_2^0,...,x_k^0)+\sum\limits_{i=1}^{k}\frac{\partial f_1}{\partial x_i^0}x_i+\\ +\frac{1}{2}\sum\limits_{i=1}^{k}\sum\limits_{j=1}^{k}\frac{\partial^2 f_1}{\partial x_i^0\partial x_j^0}x_ix_j+ o(|| X-X^0||^2),$$

где $$\lim\limits_{x\to x^0} o(|| X-X^0||^2)=0$$.

Как известно [10], вектор-градиент функции $$f_1(x_1^0,x_2^0,...,x_k^0)$$ в заданной точке $$(x_1^0,x_2^0,...,x_k^0)$$ есть

$$grad\mbox{ }f_1(x_1^0,x_2^0,...,x_k^0)= \left( \begin{array}{c} \frac{\partial f_1}{\partial x_1^0} \\ \frac{\partial f_1}{\partial x_2^0}\\ \vdots\\ \frac{\partial f_1}{\partial x_k^0} \end{array} \right) = \left( \frac{\partial f_1}{\partial x_i^0}, \frac{\partial f_1}{\partial x_i^0}, ..., \frac{\partial f_1}{\partial x_i^0} \right)^T.$$

Частные производные в (10.3) и в (10.4) рассчитываются в заданных точках и, таким образом, они представляют собой числовые значения.

Используем обозначения

$$\beta_0^0=f_1(x_1^0,x_2^0,...,x_k^0);\\ \beta_i^0=\frac{\partial f_1}{\partial x_i^0};\\ \beta_{ij}^0=\frac{\partial^2 f_1}{\partial x_i^0x_j^0};\\ \beta_{ii}^0=\frac{\partial^2 f_1}{2\partial (x_i^0)^2};\\ ......................\\$$

С учетом введенных обозначений (10.5) можно записать, что функция отклика будет определяться в следующем виде:

$$y=\beta_0^0+ \sum\limits_{i=1}^{k}\beta_i^0x_i+ \sum\limits_{1\le i<j\le k}^{}\beta_{ij}^0x_ix_j+ \sum\limits_{i=1}^{k}\beta_{ii}^0x_i^2+...+o(|| X-X^0||^2).$$

Отметим также, что

$$grad\mbox{ }f_1(X^0)=(\beta_1^0,\beta_2^0,...,\beta_k^0)^T$$

Разложение (10.6) можно рассматривать как многофакторное уравнение регрессии, линейное относительно своих коэффициентов.

Таким образом, для оценки градиента — оценки коэффициентов уравнения регрессии $$\beta_1^0,\beta_2^0,...,\beta_k^0$$ можно использовать метод наименьших квадратов, с помощью которого возможно определить все коэффициенты аппроксимации поверхности отклика и затем выделить оценки коэффициентов $$\beta_1^0,\beta_2^0,...,\beta_k^0$$, оценку градиента. Вид аппроксимации, количество слагаемых в разложении (10.6) задает исследователь по своим априорным сведениям.

2. Метод Бокса и Уильсона

Пусть функция отклика относительно кодированных переменных имеет вид

$$y=f_1(x_1,x_2,...,x_k).$$

Для функции отклика (10.7) может быть найдена оценка градиента в точке $$(x_1^0,x_2^0,...,x_k^0)$$.

Для поиска максимума функции отклика делается некоторый шаг (первый шаг) из точки $$X^0=(x_1^0,x_2^0,...,x_k^0)$$ в направлении оценки градиента $$grad\mbox{ }f_1(X^0)$$:

$$X_0^1=X^0+\alpha_0^1 grad\mbox{ }f_1(X^1_0)=\alpha_0^1 grad\mbox{ }f_1(X^0)=\alpha_0^1 \hat\beta^0$$

где:

$$\alpha_0^1 >0$$ — скаляр, параметр шага $$\Delta X_0^1=X_0^1-X^0$$ ;

$$\hat\beta^0=(\hat\beta^0_1,\hat\beta^0_2,...,\hat\beta^0_k)^T$$ — оценка градиента на предыдущем шаге;

$$X^1_0=(x_1^{10},x_2^{10},...,x_k^{10})^T$$ — новая точка факторного пространства.

Векторная запись (10.8) в координатной форме имеет вид

$$x_i^{10}=\alpha_0^1 \hat\beta_i^0,\qquadi=1,2,...,k.$$

Переход от кодированных переменных $$x_i^{10}$$ к натуральным переменным $$X_i^{10}$$ осуществляется по формуле

$$X_i^{10}=X_i^{0}+x_i^{10}S_i^{0}$$

где $$S_i^{0}=(X_{i2}^{0}-X_{i1}^{0})/2$$.

В точке кодированных переменных (факторов) $$X_{0}^{1}$$ следует произвести наблюдения функции отклика $$y_1^{10},y_2^{10},..,y_n^{10}$$ и найти ее оценку:

$$\hat y_0^1=\frac{1}{n}\sum\limits_{u=1}^{n}y_u^{10}.$$

Оценка функции отклика (10.11) есть результат измерений, которые производятся в реальных условиях с ошибками. Так как задача заключается в определении максимума функции отклика, естественно предположить, что оценка (10.11) значимо больше оценки в начальной точке, т. е. $$\hat y_0^1 >\hat y_0$$, где $$\hat y_0$$ — оценка функции отклика в точке факторного пространства $$X^0$$. Тогда можно сделать следующий шаг в направлении градиента (его оценки) $$grad\mbox{ }f_1(X^0)$$. Записывая (10.8) для $$m$$ -го шага, будем иметь точку факторного пространства, полученную в направлении оценки градиента, т. е.

$$x_0^m=\alpha_0^m\hat \beta^0,$$

где $$\alpha_0^m >\alpha_0^{m-1}>...>\alpha_0^1$$.

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

В направлении градиента производится средняя оценка функции отклика для значений $$\alpha_0^m >\alpha_0^{m-1}>...>\alpha_0^1$$ параметров шага. Как только будет найдена максимальная оценка функции отклика, то ее и принимают за оценку максимума функции отклика при движении в направлении оценки градиента [1]. Определение максимума функции отклика осуществляется простым перебором и сравнением, что, однако, приводит к увеличению затрат на проведение эксперимента.

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

Таким образом, организуется крутое восхождение к экстремуму функции отклика. В этом состоит схематичное изложение метода Бокса и Уильсона. Следует иметь в виду, что аналитический вид функции отклика, как правило, неизвестен.

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

Практическая часть

Рассмотрим процесс моделирования крутого восхождения (наискорейшего подъема) к экстремуму функции отклика. Для методической ясности проведения активного эксперимента зададим аналитический вид функции двух переменных

$$f(x,y)=3(1-x)^2 e^{(-x^2-(y+1)^2)},$$

где $$x, y$$ — независимые переменные, факторы.

Графический образ функции двух переменных (10.13) показан на рис. 10.1.

(рис 10.1) Исследуемая функция отклика

Необходимые измерения функции отклика при движении по градиенту будем проводить по аналитическому выражению (10.13).

В каждой точке факторного пространства заданную функцию отклика (10.13) будем аппроксимировать линейной моделью — линейным уравнением регрессии вида

$$f(x_1,x_2)=b_0+b_1x_1+b_2x_2,$$

где $$x_{1},\mbox{ }x_{2}$$ — факторы.

Коэффициенты $$b_{0},\mbox{ }b_{1},\mbox{ }b_{2}$$ в (10.14) определяются в постановке двухфакторного эксперимента.

Программный код моделирования крутого восхождения к экстремуму функции отклика может быть следующим:

function active;
clc,close all
syms x y
zz = 3*(1-x)^2*exp(-x^2 - (y+1)^2);
fprintf('\n\t Аналитический вид функции отклика:\n')
fprintf('\t f(x,y) = %s\n',char(zz))

% РАСЧЕТ ОПТИМАЛЬНЫХ КООРДИНАТ ФУНКЦИИ ОТКЛИКА
xy = fminsearch(@pan1,[0;0]);
fprintf('\t Оптимальные значения факторов:\n')
fprintf('\t x_optim = %g\n\t y_optim = %g\n',xy(1),xy(2))
y_max = abs(pan1(xy));
fprintf('\t Значение макcимума аналитической функции отклика:\n')
fprintf('\t f(x,y)_max = %g\n',y_max)
% ПОСТРОЕНИЕ ФУНКЦИИ ОТКЛИКА В ПРОСТРAНСТВЕ
[x1,x2] = meshgrid(-2:0.02:2,-2:0.02:1);
z =  3*(1-x1).^2.*exp(-x1.^2 - (x2+1).^2);
meshc(x1,x2,z);
grid on;
st='\bf\fontname{times new roman}\fontsize{12} Функция отклика';
title(st);
xlabel('\bf\fontname{times new roman}\fontsize{12}  x'),
ylabel('\bf\fontname{times new roman}\fontsize{12}  y'),
zlabel('\bf  f(x,y)'),
zlim([0, 4])
disp('---------------------------------------------------')
fprintf('\t ПОСТАНОВКА  ДВУХФАКТОРНОГО  ЭКСПЕРИМЕНТА\n')
fprintf('\t Уравнение линейной модели:\n')
fprintf('\t f(x1,x2) = b0 + b1*x1 + b2*x2\n')
% ------ КРУТОЕ ВОСХОЖДЕНИЯ К ЭКСТРЕМУМУ ФУНКЦИИ ОТКЛИКА ------
fprintf('\n\t 1-й ЦИКЛ РАСЧЕТА\n')
% ОСНОВНОЙ УРОВЕНЬ ФАКТОРОВ
X10 = -0.3;
X20 = -1.5;
% РАСЧЕТ ФУНКЦИИ ОТКЛИКА В ЦЕНТРЕ НАЧАЛЬНОГО ПЛАНА
f0 = abs(pan1([X10,X20]));
% ИНТЕРВАЛЫ ВАРЬИРОВАНИЯ ФАКТОРОВ
d1 = -0.2;
d2 = -0.2;
% ЗНАЧЕНИЯ ФАКТОРОВ В ТОЧКАХ НАЧАЛЬНОГО ПЛАНА
X1up = X10 + d1;
X1dw = X10 - d1;
X2up = X20 + d2;
X2dw = X20 - d2;
% ВЫВОД СООБЩЕНИЙ
fprintf('\t Координаты центра начального плана:\n\t x10 = %g\n\t x20 = %g\n',X10,X20)
fprintf('\t Интервалы варьирования:\n\t Для 1-го фактора: %g\n\t Для 2-го фактора: %g\n', d1, d2)
fprintf('\t Значения отклика в точках начального плана:\n')
y0 = [abs(pan1([X1up,X2up])); abs(pan1([X1dw,X2up]));
    abs(pan1([X1up,X2dw]));  abs(pan1([X1dw,X2dw]))];
fprintf('\t\t %g\n',y0)
% ФОРМИРОВАНИЕ МАТРИЦЫ ПЛАНИРОВАНИЯ ФАКТОРОВ
fn = rot90(ff2n(2),2);
for J = 1:length(fn(:))
    if (~fn(J)) %%% когда fn(J) == 0
        fn(J) = -1;
    
else
        continue;
    end
end
X0 = [ones(length(fn(:,1)),1),fn];
%----------------------------------------------
fprintf('\tМатрица планирования двухфакторного эксперимента:\n')
fprintf('\t x0\t\t\tx1\t\t  x2\t\t y\n')
disp([X0,y0])
%----------------------------------------------
% РАСЧЕТ КОЭФФИЦИЕНТОВ ЛИНЕЙНОЙ МОДЕЛИ
b0 = regress(y0,X0);
fprintf('\t Коэффициенты линейной модели:\n')
for J = 1:length(b0) 
    fprintf('\t\t b%d = %g\n',J-1,b0(J))
end
%----------------------------------------------
fprintf('\t Значения градиента:\n')
for J = 2:length(b0)
    fprintf('\t g%d = %g\n',J-1,b0(J))
end
%----------------------------------------------
% ОРГАНИЗАЦИЯ ДВИЖЕНИЯ ПО ГРАДИЕНТУ
k = 0;
for J = 0.2:0.2:4
    k = k + 1;
aa(k) = J;
x11 = J*b0(2);
x22 = J*b0(3);
X1(k) = X10 + x11*d1;
X2(k) = X20 + x22*d2;
f1(k) = abs(pan1([X1(k),X2(k)]));
end
%-----------------------------------------------
k1 = find(f1 == max(f1));
% КООРДИНАТЫ ЦЕНТРА НОВОГО ПЛАНА
x10 = X1(k1);
x20 = X2(k1);
%-----------------------------------------------
fprintf('\t Число шагов поиска максимума функции отклика по направлению градиента: k = %d\n',k1)
fprintf('\t Параметр шага в направлении оценки градиента: a = %g\n',aa(k1))
fprintf('\t Расчетное значение функции отклика на градиенте: y1 = %g\n',max(f1))
fprintf('\t Координаты центра нового плана:\n\t x11 = %g\n\t x21 = %g\n',X1(k1),X2(k1))
disp('-------------------------------------------------')

fprintf('\t 2-й ЦИКЛ РАСЧЕТА\n')
% КООРДИНАТЫ ФАКТОРОВ НОВОГО ПЛАНА

X11up = x10 + d1;
X11dw = x10 - d1;
X21up = x20 + d2;
X21dw = x20 - d2;
%----------------------------------------------
fprintf('\t Новые интервалы варьирования:\n\t Для 1-го фактора: %g\n\t Для 2-го фактора: %g\n', d1, d2)
fprintf('\t Значения отклика в точках нового плана:\n')
y01 = [abs(pan1([X11up,X21up])); abs(pan1([X11dw,X21up]));
       abs(pan1([X11up,X21dw])); abs(pan1([X11dw,X21dw]))];
fprintf('\t\t %g\n',y01)
%----------------------------------------------
% ФОРМИРОВАНИЕ МАТРИЦЫ ПЛАНИРОВАНИЯ
fn = rot90(ff2n(2),2);
for J = 1:length(fn(:))
    if (~fn(J)) %% fn(J) == 0
        fn(J) = -1;
    else
        continue;
    end
end
X0 = [ones(length(fn(:,1)),1),fn];
%----------------------------------------------
fprintf('\tМатрица планирования двухфакторного эксперимента:\n')
fprintf('\t x0\t\t\tx1\t\t  x2\t\t y\n')
disp([X0,y01])

% РАСЧЕТ КОЭФФИЦИЕНТОВ ЛИНЕЙНОЙ МОДЕЛИ
b01 = regress(y01,X0);
fprintf('\t Коэффициенты линейной модели:\n')
for J = 1:length(b01) 
    fprintf('\t\t b%d = %g\n',J-1,b01(J))
end
%-----------------------------------------------
fprintf('\t Значения градиента:\n')
for J = 2:length(b01)
    fprintf('\t g%d = %g\n',J-1,b01(J))
end
%-----------------------------------------------
% ОРГАНИЗАЦИЯ ДВИЖЕНИЯ ПО ГРАДИЕНТУ
k = 0;
for J = 0.1:0.1:4
    k = k + 1;
aa(k) = J;
x11 = J*b01(2);
x22 = J*b01(3);
X11(k) = x10 + x11*d1;
X21(k) = x20 + x22*d2;
f11(k) = abs(pan1([X11(k),X21(k)]));
end
%----------------------------------------------

k1 = find(f11 == max(f11));
% КООРДИНАТЫ ЦЕНТРА НОВОГО ПЛАНА
x110 = X11(k1);
x220 = X21(k1);
%----------------------------------------------
fprintf('\t Число шагов поиска максимума функции отклика по направлению градиента: k = %d\n',k1)
fprintf('\t Параметр шага в направлении оценки градиента: a = %g\n',aa(k1))
fprintf('\t Расчетное значение функции отклика на градиенте: y2 = %g\n',max(f11))
fprintf('\t Координаты центра нового плана:\n\t x11 = %g\n\t x21 = %g\n',X11(k1),X21(k1))
disp('-------------------------------------------------')

fprintf('\t 3-й ЦИКЛ РАСЧЕТА\n')
% КООРДИНАТЫ ФАКТОРОВ НОВОГО ПЛАНА
X11up = x110 + d1;
X11dw = x110 - d1;
X21up = x220 + d2;
X21dw = x220 - d2;
%----------------------------------------------
fprintf('\t Новые интервалы варьирования:\n\t Для 1-го фактора: %g\n\t Для 2-го фактора: %g\n', d1, d2)
fprintf('\t Значения отклика в точках нового плана:\n')
y03 = [abs(pan1([X11up,X21up])); abs(pan1([X11dw,X21up]));
       abs(pan1([X11up,X21dw])); abs(pan1([X11dw,X21dw]))];
fprintf('\t\t %g\n',y01)
%----------------------------------------------
% ФОРМИРОВАНИЕ МАТРИЦЫ ПЛАНИРОВАНИЯ
fn = rot90(ff2n(2),2);
for J = 1:length(fn(:))
    if (~fn(J)) %% fn(J) == 0
        fn(J) = -1;
    else
        continue;
    end
end
X0 = [ones(length(fn(:,1)),1),fn];
%---------------------------------------------
fprintf('\tМатрица планирования двухфакторного эксперимента:\n')
fprintf('\t x0\t\t\tx1\t\t  x2\t\t y\n')
disp([X0,y03])
%---------------------------------------------
% РАСЧЕТ КОЭФФИЦИЕНТОВ ЛИНЕЙНОЙ МОДЕЛИ
b03 = regress(y03,X0);
fprintf('\t Коэффициенты линейной модели:\n')
for J = 1:length(b03) 
    fprintf('\t\t b%d = %g\n',J-1,b03(J))
end
%---------------------------------------------

fprintf('\t Значения градиента:\n')
for J = 2:length(b03)
    fprintf('\t g%d = %g\n',J-1,b03(J))
end
%----------------------------------------------
% ОРГАНИЗАЦИЯ ДВИЖЕНИЯ ПО ГРАДИЕНТУ
k = 0;
for J = 0.2:0.2:2
    k = k + 1;
aa(k) = J;
x11 = J*b03(2);
x22 = J*b03(3);
X11(k) = x110 + x11*d1;
X21(k) = x220 + x22*d2;
f33(k) = abs(pan1([X11(k),X21(k)]));
end
%-----------------------------------------------
k3 = find(f33 == max(f33));
% КООРДИНАТЫ ЦЕНТРА НОВОГО ПЛАНА
x110 = X11(k3);
x220 = X21(k3);
%------------------------------------------------
fprintf('\t Число шагов поиска максимума функции отклика по направлению градиента: k = %d\n',k3)
fprintf('\t Параметр шага в направлении оценки градиента: a = %g\n',aa(k3))
fprintf('\t Расчетное значение функции отклика на градиенте: y3 = %g\n',max(f33))
fprintf('\t Координаты центра нового плана:\n\t x11 = %g\n\t x21 = %g\n',X11(k3),X21(k3))
disp('-------------------------------------')

fprintf('\t 4-й ЦИКЛ РАСЧЕТА\n')
% ИНТЕРВАЛЫ ВАРЬИРОВАНИЯ ФАКТОРОВ
d41 = -0.1; d42 = -0.1;
% КООРДИНАТЫ ФАКТОРОВ НОВОГО ПЛАНА
X11up = x110 + d41;
X11dw = x110 - d41;
X21up = x220 + d42;
X21dw = x220 - d42;
%-------------------------------------------

fprintf('\t Новые интервалы варьирования\n\t Для 1-го фактора: %g\n\t Для 2-го фактора: %g\n', d41, d42)
fprintf('\t Значения отклика в точках нового плана:\n')
y04 = [abs(pan1([X11up,X21up])); abs(pan1([X11dw,X21up]));
       abs(pan1([X11up,X21dw])); abs(pan1([X11dw,X21dw]))];
fprintf('\t\t %g\n', y04)

% ФОРМИРОВАНИЕ МАТРИЦЫ ПЛАНИРОВАНИЯ
fn = rot90(ff2n(2), 2);

for J = 1:length(fn(:))
    if (~fn(J))   %% когда fn(J) == 0
        fn(J) = -1;
    else
        continue;
    end
end
X0 = [ones(length(fn(:,1)),1), fn];
%---------------------------------------------
fprintf('\tМатрица планирования двухфакторного эксперимента:\n')
fprintf('\t x0\t\t\tx1\t\t  x2\t\t y\n')
disp([X0, y04])
%----------------------------------------------
% РАСЧЕТ КОЭФФИЦИЕНТОВ ЛИНЕЙНОЙ МОДЕЛИ
b04 = regress(y04,X0);
fprintf('\t Коэффициенты линейной модели:\n')
for J = 1:length(b04) 
    fprintf('\t\t b%d = %g\n', J-1, b04(J))
end
fprintf('\t Значения градиента:\n')
for J = 2:length(b04)
    fprintf('\t g%d = %g\n', J-1, b04(J))
end
%----------------------------------------------
% ОРГАНИЗАЦИЯ ДВИЖЕНИЯ ПО ГРАДИЕНТУ
k = 0;
for J = 5:0.5:10
    k = k + 1;
aa(k) = J;
x11 = J*b04(2);
x22 = J*b04(3);
X11(k) = x110 + x11*d41;
X21(k) = x220 + x22*d42;
f44(k) = abs(pan1([X11(k), X21(k)]));
end
%-------------------------------------------------
k4 = find(f44 == max(f44));
% КООРДИНАТЫ ЦЕНТРА НОВОГО ПЛАНА
x110 = X11(k4);
x220 = X21(k4);
%-------------------------------------------------

fprintf('\t Число шагов поиска максимума функции отклика по направлению градиента: k = %d\n', k4)
fprintf('\t Параметр шага в направлении оценки градиента: a = %g\n',aa(k4))
fprintf('\t Расчетное значение функции отклика на градиенте: y4 = %g\n',max(f44))
fprintf('\t Координаты центра нового плана:\n\t x11 = %g\n\t x21 = %g\n',X11(k4),X21(k4))
disp('-----------------------------------------------')
 
function f = pan1(w)
%%% М-файл описания заданной функции
%%% двух переменных
x = w(1);
y = w(2);
f = 3*(1-x)^2*exp(-x^2-(y+1)^2);

Результат выполнения программы

Аналитический вид функции отклика:
	 f(x,y) = 3*(1-x)^2*exp(-x^2-(y+1)^2)
	 Оптимальные значения факторов:
	 x_optim = -0.618038
	 y_optim = -0.999985
	 Значение макcимума аналитической функции отклика:
	 f(x,y)_max = 5.36057
--------------------------------------------------------------
	 ПОСТАНОВКА  ДВУХФАКТОРНОГО  ЭКСПЕРИМЕНТА
	 Уравнение линейной модели:
	 f(x1,x2) = b0 + b1*x1 + b2*x2
	 1-й ЦИКЛ РАСЧЕТА
	 Координаты центра начального плана:
	 x10 = -0.3
	 x20 = -1.5
	 Интервалы варьирования:
	 Для 1-го фактора: -0.2
	 Для 2-го фактора: -0.2
	 Значения отклика в точках начального плана:
		 3.22052
		 2.20171
		 4.80445
		 3.28456
	Матрица планирования двухфакторного эксперимента:
      x0		x1	     x2		 y
    1.0000    1.0000    1.0000    3.2205
    1.0000   -1.0000    1.0000    2.2017
    1.0000    1.0000   -1.0000    4.8044
    1.0000   -1.0000   -1.0000    3.2846
	 Коэффициенты линейной модели:
		 b0 = 3.37781
		 b1 = 0.634676
		 b2 = -0.666696
	 Значения градиента:
	 g1 = 0.634676
	 g2 = -0.666696
Число шагов поиска максимума функции отклика по направлению градиента: k = 15
	 Параметр шага в направлении оценки градиента: a = 3
	 Расчетное значение функции отклика на градиенте: y1 = 5.27863
	 Координаты центра нового плана:
	 x11 = -0.680805
	 x21 = -1.09998
---------------------------------------------------------------
	 2-й ЦИКЛ РАСЧЕТА
	 Новые интервалы варьирования:
	 Для 1-го фактора: -0.2
	 Для 2-го фактора: -0.2
	 Значения отклика в точках нового плана:
		 4.46471
		 4.77131
		 4.8365
		 5.16863
	Матрица планирования двухфакторного эксперимента:
     x0	       x1	     x2		 y
    1.0000    1.0000    1.0000    4.4647
    1.0000   -1.0000    1.0000    4.7713
    1.0000    1.0000   -1.0000    4.8365
    1.0000   -1.0000   -1.0000    5.1686
	 
      Коэффициенты линейной модели:
		 b0 = 4.81029
		 b1 = -0.159682
		 b2 = -0.192275
	 
     Значения градиента:
	 g1 = -0.159682
	 g2 = -0.192275
Число шагов поиска максимума функции отклика по направлению градиента: k = 23
	 Параметр шага в направлении оценки градиента: a = 2.3
	 Расчетное значение функции отклика на градиенте: y2 = 5.35901
	 Координаты центра нового плана:
	 x11 = -0.607352
	 x21 = -1.01154
-----------------------------------------------------------------
	 3-й ЦИКЛ РАСЧЕТА
	 Новые интервалы варьирования:
	 Для 1-го фактора: -0.2
	 Для 2-го фактора: -0.2
	 Значения отклика в точках нового плана:
		 4.46471
		 4.77131
		 4.8365
		 5.16863
	Матрица планирования двухфакторного эксперимента:
      x0	      x1         x2		 y
    1.0000    1.0000    1.0000    4.8831
    1.0000   -1.0000    1.0000    4.8131
    1.0000    1.0000   -1.0000    4.9283
    1.0000   -1.0000   -1.0000    4.8578
	 Коэффициенты линейной модели:
		 b0 = 4.87058
		 b1 = 0.0351275
		 b2 = -0.0224743
	 

     Значения градиента:
	 g1 = 0.0351275
	 g2 = -0.0224743
Число шагов поиска максимума функции отклика по направлению градиента: k = 9
	 Параметр шага в направлении оценки градиента: a = 1.8
	 Расчетное значение функции отклика на градиенте: y3 = 5.36048
	 Координаты центра нового плана:
	 x11 = -0.619997
	 x21 = -1.00345
-----------------------------------------------------------------
	 4-й ЦИКЛ РАСЧЕТА
	 Новые интервалы варьирования:
	 Для 1-го фактора: -0.1
	 Для 2-го фактора: -0.1
	 
      Значения отклика в точках нового плана:
		 5.22869
		 5.23272
		 5.2359
		 5.23994
	
   Матрица планирования двухфакторного эксперимента:
      x0 	      x1         x2		 y
    1.0000    1.0000    1.0000    5.2287
    1.0000   -1.0000    1.0000    5.2327
    1.0000    1.0000   -1.0000    5.2359
    1.0000   -1.0000   -1.0000    5.2399
	 Коэффициенты линейной модели:
		 b0 = 5.23431
		 b1 = -0.002017
		 b2 = -0.00360653
	 Значения градиента:
	 g1 = -0.002017
	 g2 = -0.00360653
Число шагов поиска максимума функции отклика по направлению градиента: k = 10
	 Параметр шага в направлении оценки градиента: a = 9.5
	 Расчетное значение функции отклика на градиенте: y4 = 5.36057
	 Координаты центра нового плана:
	 x11 = -0.618081
	 x21 = -1.00002
-----------------------------------------------------------------

Задание

На поверхности функции отклика вычертите линии движения к экстремуму.

Выведите значения функции отклика в точках оптимального плана, когда в центре его достигается оценка максимума функции отклика.

Примените функцию $$ff2n$$ (см. $$help\mbox{ }ff2n$$ ) без перестановки столбцов, которая в программе выполнялпсь с помощью $$rot90$$. Проанализируйте полученные результаты.

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

Напишите программу поиска максимума функции отклика с автоматическим выбором количества циклов, значений шага по градиенту, с заданием в диалоговых окнах $$inputdlg$$ интервалов варьирования и центра начального плана.

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

Что такое функция отклика?

В чем различие между пассивным и активным экспериментами?

Какие задачи решает активный эксперимент?

Что называется градиентом функции отклика?

В чем заключается методика поиска экстремума функции отклика?

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

Как осуществляется построение факторного эксперимента в окрестности заданной точки факторного пространства?

На каком этапе поиска экстремума функции отклика применяется метод наименьших квадратов?