1-ая лабораторная - 16-ый вариант - решение, код задачи по моделированию систем - ФКН ВГУ

ПРИМЕЧАНИЕ: вместо мат ожидания проведено исследование дисперсии, а потому условие задачи переписано таким образом:

16. Провести стратегическое и тактическое планирование модельного эксперимента. Выходной реакцией системы является случайная величина, распределенная по закону Эрланга. Факторами являются параметры: b (4; 6), c (3; 6). Оценить показатель эффективности системы – дисперсию реакции. Доверительный интервал dm = 0,2 с уровнем значимости ? = 0,05.

Код (с комментариями):

clear all;
clear functions;

nf=2;
minf=[4 3];%minf=[1 0.5];
maxf=[6 6];%maxf=[5 1];
%формирование дробного двухуровневого плана эксперимента
% Формирование плана фактически также используется для СТРАТЕГИЧЕСКОГО
% планирования
%для учета взаимодействий
fracfact('a b ab' );
N=2^nf;
fracplan=ans
fictfact=ones(N,1); % фиктивный фактор
X = [fictfact ans]'
fraceks = zeros(N,nf);
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

%тактическое планирование эксперимента
%задание доверительного интервала и уровня значимости
d_sigma=0.2;
alpha=0.05;
%определение 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);
% вложенный цикл статистических испытаний (ТАКТИЧЕСКОЕ ПЛАНИРОВАНИЕ)
    for k=1:NE,
        % тактическое планирование необходимо чтобы проверить устойчивость
        % системы к всевозможным случайным событиям (напр. событиям внешней среды)
        %имитация функционирования системы
        u(k)=systemqv(a,b);% пользовательская функция реализующая датчик закона распределения
%она определяется так (в моём случае это Эрлагановское распределение):
        % function u = systemqv(b,c)
        %      product=1
        %      for i=1:c
        %            product=product*randn
        %      end
        %      u=(-b)*log(product);
        %  return;
    end; % конец цикла тактического планирования
    % далее мы  - после формирования массива можем посчитать либо
    % мат.ожидание либо дисперсию
    %оценка параметров (реакции) по выборке наблюдений
    % ДАЛЕЕ для построения поверхности мы используем фиктивный фактор  -
    % чтобы свернуть формулу 2.1 в формулу 2.2 ()
        mx=mean(u);
        DX=std(u)^2;
    Y(j)=DX; % сюда пишем полученный результат для мат ожидания (дисперсии)
    % здесь мы заканчиваем проводить первый эксперимент - данные действия
    % надо повторить для каждого эксперимента
    %формирование и отображение гистограммы с 12-ю интервалами
    %figure;
    %hist(u,12);
end; % конец цикла стратегического планирования


% закончили проводить эксперимент и 
% далее занялись построением поверхностей - экспериментальная 
% строится на основе эксперимента (что логично)
%определение коэффициентов регрессии
C=X*X';
b_=inv(C)*X*Y'
%формирование зависимости реакции системы на множестве
%значений факторов
A=minf(1):0.1:maxf(1); % 0.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);% формула 2.2 (Стратегическое планирование.pdf)
    end;
end;
for i=1:N1,
    for j=1:N2,
        %реальная поверхность реакции
            %Yo(j,i)=A(i)*B(j);
            Yo(j,i)=(A(i)^2)*B(j);
            %Yo(j,i)=(A(i)^2)*exp(B(j)^2)*(exp(B(j)^2)-1);
    end;
end;

% далее - РИСУЕМ ГРАФИКИ (теоретический и экспериментальный)
% отображение зависимостей в трехмерной графике 
[x,y]=meshgrid(A,B);
figure;
% рисуем график полученный в ходе эксперимента
subplot(1,2,1),plot3(x,y,Yc),% 1 ячейка по горизонтали - 2 по вертикали
xlabel('fact a'),
ylabel('fact b'),
zlabel('Yc'),
title('System output'),
grid on,

% далее рисуем "реальный" график
subplot(1,2,2),plot3(x,y,Yo),
xlabel('fact a'),
ylabel('fact b'),
zlabel('Yo'),
title('System output'),
grid on;