Код и комментарии к первой задачи (пример) по Моделированию систем ФКН ВГУ 4 курс 2013

%Стратегическое и тактическое планирование модельного эксперимента
clear all;
%задание количества факторов и диапазонов значений факторов
nf=3; % задаётся количество факторов, которое мы будем рассматривать
minf=[-5 -10 0]; % первый фактор от -5 до 5  - второй (-10,10)
maxf=[5 10 20]; %третий (0 - 20)
%задание количества уровней каждого фактора
level=[3 3 2]; % первый фактор 
% уровень - показывает сколько различных значений может принимать фактор в
% заданном диапазоне - но не одновременно, а вообще


% мы используем матлаб для данной лабы так как здесь куча готовых функций
%формирование полного плана эксперимента 
fullfact(level);% даёт набор некоторыз значений, 
% которые по сути являются факторным планом (все возможные значения для всех факторов)
%ans - хранит результат предыдущей операции
fulplan = ans
% запустим весь код выше (до fulplan = ans)
% - первый и второй столбец изменяются 123-123, а
% третий - только 2   =
%      2     2     1
%      3     2     1
%      1     3     1
%      2     3     1
%      3     3     1
%      1     1     2
%      2     1     2
%      3     1     2
%      1     2     2
%      2     2     2
%      3     2     2
%      1     3     2
%      2     3     2
%      3     3     2
%      значение - 1 или 2 (но это лишь абстрактыные значения)
% Далее в цикле мы заменяем абстрактные значения на реальные
N=3*3*2; %число возможных перестановок (1)
for i=1:nf,
    for j=1:N,
fuleks(j,i)=minf(i)+(fulplan(j,i)-1)*(maxf(i)-minf(i))/(level(i)-1);
end;
end;
fuleks
% запустим код до предыдущей строки от строки (1)(предыдущие 6 строк)
% и получим массив:
%      0     0     0
%      5     0     0
%     -5    10     0
%      0    10     0
%      5    10     0
%     -5   -10    20
%      0   -10    20
%      5   -10    20
%     -5     0    20
%      0     0    20
%      5     0    20
%     -5    10    20
%      0    10    20
%      5    10    20
% Это уже реальные значения из диапазона

% Можно при стратегическом планировани рассматривать только граничные точки
% - для формирования такого плана жксперимента используется 
%формирование полного двухуровневого плана эксперимента использует ff2n -
%здесь(далее) у нас фиксированно рассматриваются два уровня эксперимента
ff2n(nf);
N=2^nf;
ff2nplan=ans
for i=1:nf,
    for j=1:N,
fuleks2n(j,i)=minf(i)+ff2nplan(j,i)*(maxf(i)-minf(i));
end;
end;
fuleks2n
% запускаем предыдущие строки
% fuleks2n =
% 
%     -5   -10     0
%     -5   -10    20
%     -5    10     0
%     -5    10    20
%      5   -10     0
%      5   -10    20
%      5    10     0
%      5    10    20
%формирование дробного двухуровневого плана эксперимента
% Далее мы рассматрием не просто значение каждого фактора, но и их
% взаимодействия
N=2^nf;
fracfact('a b c ab bc ac abc' ); % список возможных взимодействий
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 % ДО СЮДА РАЗБИРАЕМ ДЛЯ ПЕРВОЙ ЗАДАЧИ
% запустив предыдущее получим
% fracplan =
% Итак сначала мы задаём =  a b c ab bc ac   abc (столбец для каждого)
%     -1    -1    -1     1     1     1    -1
%     -1    -1     1     1    -1    -1     1
%     -1     1    -1    -1    -1     1     1
%     -1     1     1    -1     1    -1    -1
%      1    -1    -1    -1     1    -1     1
%      1    -1     1    -1    -1     1    -1
%      1     1    -1     1    -1    -1    -1
%      1     1     1     1     1     1     1
% 
% 
% fraceks =
% 
%     -5   -10     0
%     -5   -10    20
%     -5    10     0
%     -5    10    20
%      5   -10     0
%      5   -10    20
%      5    10     0
%      5    10    20
% до этого места всё выше было закомментированно
clear all;
nf=2;
minf=[1 0.5];
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.1;
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)=systemeqv(a,b);% пользовательская функция реализующая датчик закона распределения
%она определяется так:
        %         function u=systemeqv(a,b);
%           %датчик логнормального распределения
%           u=a*exp(b*randn);
%           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)^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;
% fulplan =
% 
%      1     1     1
%      2     1     1
%      3     1     1
%      1     2     1
%      2     2     1
%      3     2     1
%      1     3     1
%      2     3     1
%      3     3     1
%      1     1     2
%      2     1     2
%      3     1     2
%      1     2     2
%      2     2     2
%      3     2     2
%      1     3     2
%      2     3     2
%      3     3     2
% 
% 
% fuleks =
% 
%     -5   -10     0
%      0   -10     0
%      5   -10     0
%     -5     0     0
%      0     0     0
%      5     0     0
%     -5    10     0
%      0    10     0
%      5    10     0
%     -5   -10    20
%      0   -10    20
%      5   -10    20
%     -5     0    20
%      0     0    20
%      5     0    20
%     -5    10    20
%      0    10    20
%      5    10    20
% 
% 
% ff2nplan =
% 
%      0     0     0
%      0     0     1
%      0     1     0
%      0     1     1
%      1     0     0
%      1     0     1
%      1     1     0
%      1     1     1
% 
% 
% fuleks2n =
% 
%     -5   -10     0
%     -5   -10    20
%     -5    10     0
%     -5    10    20
%      5   -10     0
%      5   -10    20
%      5    10     0
%      5    10    20
% 
% 
% fracplan =
% 
%     -1    -1    -1     1     1     1    -1
%     -1    -1     1     1    -1    -1     1
%     -1     1    -1    -1    -1     1     1
%     -1     1     1    -1     1    -1    -1
%      1    -1    -1    -1     1    -1     1
%      1    -1     1    -1    -1     1    -1
%      1     1    -1     1    -1    -1    -1
%      1     1     1     1     1     1     1
% 
% 
% fraceks =
% 
%     -5   -10     0
%     -5   -10    20
%     -5    10     0
%     -5    10    20
%      5   -10     0
%      5   -10    20
%      5    10     0
%      5    10    20
% 
% 
% fracplan =
% 
%     -1    -1     1
%     -1     1    -1
%      1    -1    -1
%      1     1     1
% 
% 
% X =
% 
%      1     1     1     1
%     -1    -1     1     1
%     -1     1    -1     1
%      1    -1    -1     1
% 
% 
% fraceks =
% 
%     1.0000    0.5000
%     1.0000    1.0000
%     5.0000    0.5000
%     5.0000    1.0000
% 
% 
% NE =
% 
%    769
% 
% 
% b_ =
% 
%    26.6395
%    24.8441
%    21.6868
%    20.2838
% 
%function u=systemeqv(a,b);
%логнормальное распределение с параметрами формы a,b
%u=a*exp(b*randn);


 function u = systemqv(b,c)
      product=1
      for i=1:c
            product=product*randn
      end
      u=(-b)*log(product);
  return;
vedro-compota's picture

это к нашей задаче или вообще ?)

_____________
матфак вгу и остальная классика =)

к нашей задаче,16 и 10 вариант!

vedro-compota's picture

ок) спасибо) я сначала подумал что мы здесь о информационной безопасности - пора уже отдельный раздел форума сделать))

_____________
матфак вгу и остальная классика =)