%очистка рабочего пространства
clear all;
%задание количества факторов
nf=2;
%задание диапазонов значений факторов (минимальное и максимальное значение)
minf=[-2 0.1];
maxf=[0.6 0.6];
%формирование дробного двухуровневого плана эксперимента
%для учета взаимодействий (три столбца)
fracfact('a b ab' );
%объем обучающей выборки
N=2^nf;
fracplan=ans
%добавление фиктивного фактора чтоб была матрица 4*4
fictfact=ones(N,1);
X=[fictfact ans]'
%заполнение матрицы значениями факторов
fraceks=zeros(N,nf);
fraceks
%факторы:испытания
for i=1:nf,
    for j=1:N,
        fraceks(j,i)=minf(i)+(fracplan(j,i)+1)*(maxf(i)-minf(i))/2;
    end;
end;
fraceks
%тактическое планирование эксперимента
%задание доверительного интервала
dm=0.3;
%задание уровня значимости
alpha=0.15;
%определение t-критического
%t-критическое (из формулы, которую мы выписывали)
tkr_alpha=norminv(1-alpha/2);
%определение требуемого числа испытаний
%NE=round(1+2*tkr_alpha^2/d_sigma^2)
%цикл по совокупности экспериментов стратегического плана
for j=1:N
    a=fraceks(j,1);
    b=fraceks(j,2);
    u=[];
    %передача факторов в нашу функцию
    u(1)=systemeqv(a,b);
    NE=2;
    %счетчик
    i=2;
    %цикл повторяется пока NE>=i
    %NE уменьшается т.к. при увеличении количества точек,
    %дисперсия будет уменьшаться
    while NE>=i
        u(i)=systemeqv(a,b);
        NE=round((std(u)^2*tkr_alpha^2)/(dm^2));
        i=i+1;
    end;
    %4 крайние точки
    Y(j)=mean(u);
end;

%определение коэффициентов регрессии
C=X*X';
b_=inv(C)*X*Y'
%формирование зависимости реакции системы на множестве
%значений факторов
%минимальные шаги между двумя линиями
%0.1 - это шаг
A=minf(1):0.1:maxf(1);
B=minf(2):0.1:maxf(2);
%заполнение матриц числами
[k N1]=size(A);
[k N2]=size(B);
for i=1:N1,
    for j=1:N2,
        an(i)=2*(A(i)-minf(1))/(maxf(1)-minf(1))-1;
        bn(j)=2*(B(j)-minf(2))/(maxf(2)-minf(2))-1;
        %экспериментальная поверхность реакции
        Yc(j,i)=b_(1)+an(i)*b_(2)+bn(j)*b_(3)+an(i)*bn(j)*b_(4);
end;
end;
for i=1:N1,
    for j=1:N2,
        %реальная поверхность реакции
        Yo(j,i)=A(i)*exp((B(j)^2)/2);
end;
end;
%отображение графиков с помощью соответствующих функций
[x,y]=meshgrid(A,B);
figure;
%график полученный в ходе эксперимента
subplot(1,2,1),surfl(x,y,Yc), %создание трехмерного графика  с бликами света
xlabel('fact a'),
ylabel('fact b'),
zlabel('Yc'),
title('Experimental surface output'),
grid on,
%реальный гафик
subplot(1,2,2),surfl(x,y,Yo), %создание трехмерного графика  с бликами света
xlabel('fact a'),
ylabel('fact b'),
zlabel('Yo'),
title('Real surface output'),
grid on;