Моделирование систем. Моделирование систем 2е издание, исправленное
 Скачать 5.72 Mb.
|
7. В каком случае влияние качественного фактора на выходную переменную отсутствует? 8. Как рассчитывается несмещенная оценка выборочной дисперсии? Планирование и обработка результатов пассивного эксперимента Цель работы: практически изучить применение линейных регрессионных моделей эксперимента с помощью компьютерного моделирования для случая, когда объект исследования по техническим, технологическим или экономическим соображениям не допускает преднамеренного варьирования входных переменных в необходимом диапазоне, и о виде математической модели и ее параметров делается заключение по результатам наблюдений входных и выходных переменных в режиме нормального функционирования исследуемого объекта или системы. Среда программирования — MATLAB. Теоретическая часть Основная задача пассивного эксперимента — по результатам наблюдений сделать некоторые выводы о параметрах математической модели эксперимента [1] . При этом вид ее предполагается известным, а параметры — неизвестными. Далее будет рассматриваться класс линейных регрессионных моделей эксперимента. В общем случае объект исследования представляется в виде схемы, представленной на рис. 9.1 , где — входные величины или факторы; — -я выходная величина ( ); — случайные неконтролируемые возмущения [12] Под моделью объекта по -му каналу понимают функцию (9.1) Так как имеются случайные неконтролируемые возмущения, изменение функции (9.1) носит случайный характер, а потому для получения математического описания (9.1) применяются методы регрессионного анализа на основе статистических данных, накопленных в результате проведения эксперимента. Рис. 9.1. Схема объекта исследования Применение метода пассивного эксперимента может быть успешным, если при его проведении соблюдаются необходимые условия, к которым относятся такие, как правильное определение времени регистрации данных, обеспечение независимости соседних измерений и входных переменных друг от друга, достаточный с точки зрения математической статистики объем экспериментальных данных [12] Если функция не имеет бесконечных разрывов, то ее можно разложить в степенной ряд Тейлора: (9.2) где , , , — постоянные коэффициенты уравнения, оценки которых необходимо определить в результате постановки и проведения пассивного эксперимента; — число наиболее существенных входных величин, полученных в результате отсеивающего эксперимента. Пространство, в котором строится поверхность отклика – реакция выходной величины, называется факторным пространством . Для применения методов регрессионного анализа требуется выполнение ряда предпосылок [12] , а именно: 1. результаты наблюдений выходной величины в точках факторного пространства представляют собой независимые случайные величины, распределенные по нормальному закону, а процесс изменения выходной величины должен быть стационарным во времени; 2. дисперсии этих случайных величин должны быть равны друг другу (выборочные оценки дисперсий однородны); 3. все значения входных величин должны измеряться с пренебрежимо малой ошибкой по сравнению с ошибкой измерения выходной величины; 4. входные величины не должны коррелировать между собой; 5. все соседние измерения по каждой -й входной величине должны быть независимы. Число коэффициентов уравнения (9.2) определяет объем эксперимента. Поэтому выбирают такой полином, который содержит как можно меньше коэффициентов, но удовлетворяет требованию простоты и адекватности, под которой понимается способность модели предсказывать результаты эксперимента в некоторой области с требуемой точностью. Часто на предварительной стадии исследования объекта выбирают полином первой степени (включая первую сумму в уравнении (9.2), т. е. в разложении Тейлора), предполагая, что параметры объекта лежат в области, в которой расположен экстремум исследуемой функции (выходная величина системы – отклик), и поэтому объект описывается линейной моделью. Если же эта линейная модель оказывается неадекватной, то в нее включают члены парного взаимодействия (включая вторую сумму в уравнении (9.2), т. е. в разложении Тейлора), а при необходимости увеличивают степень полинома до тех пор, пока модель не окажется адекватной [12] В результате регрессионного анализа результатов пассивного эксперимента находят оценки коэффициентов уравнения регрессии 1. План эксперимента Рассмотрим эксперимент, в котором проводится измерений зависимой переменной в некоторых точках факторного пространства [1] . Обозначим через наблюдаемое значение зависимой переменной в -м опыте в точке Здесь — значение переменной в -м опыте. Определение 1. Набор точек называется планом эксперимента. Точки при этом не обязательно должны быть различными. Матрица вида (9.3) называется матрицей плана эксперимента. Если обозначить различные точки плана через , то совокупность таких точек называется спектром плана. Определение 2. Нормированным планом называют совокупность величин ; , где В пассивном эксперименте задача построения плана не рассматривается. Матрица плана (9.3) предполагается известной (заданной) или является предопределенной условиями проведения эксперимента [1] . Задача исследователя в пассивном эксперименте состоит в выполнении наблюдений над выходной (зависимой) переменной в точках, определяемых матрицей плана, и последующем анализе их результатов. Примечание. В случае, когда регрессионная модель строится непосредственно по измерениям без предварительного усреднения, для того чтобы произвести оценку свободного члена, в модели необходимо представить матрицу плана в виде (9.4) 2. Одномерная регрессионная модель эксперимента В общем случае вид функции отклика неизвестен. Будем предполагать, что функция отклика является одномерной и представима в виде (9.5) где: есть -мерный вектор неизвестных параметров; — известные функции, которые называются еще базисными функциями. Под регрессионной моделью эксперимента будем понимать линейную по параметрам функцию отклика (9.5). В качестве примера рассмотрим функцию отклика Построим матрицу планирования по заданным измерениям входной переменной . В соответствии с (9.4) будем иметь Еще один пример. Пусть функция отклика имеет вид Значения независимых переменных равны следующим значениям: Тогда матрица планирования будет иметь такой вид, что в верхней строчке ее будут указаны аргументы функции отклика где — фиктивная переменная. 3. Оценивание параметров одномерной функции отклика В случае, когда одномерная функция отклика является линейной относительно неизвестных параметров, т. е. в произвольной точке по результатам наблюдений (измерений) в точках , где , задаваемых матрицей плана (планирования) , наилучшая оценка в смысле метода наименьших квадратов равна (9.6) если ранг матрицы планирования равен числу неизвестных параметров, т. е. , где — число неизвестных параметров [1] Оценка параметров по формуле (9.6) называется оценкой по наблюдениям полного ранга, т. е. когда . В этом случае формула (9.6) вытекает из так называемого нормального уравнения (9.7) где — информационная матрица. В методе наименьших квадратов минимизируется функционал (скалярная величина) (9.8) где — вектор неизвестных параметров, подлежащих оценке. Оценка параметров по формуле (9.6) доставляет минимум функционалу (9.8). При подстановке (9.6) в (9.8) получающаяся величина называется остаточной суммой квадратов — RSS (Residual Sum of Squares). 4. Оценивание параметров многомерной функции отклика В случае, когда имеются несколько выходных переменных — функций отклика, для каждой функции отклика можно произвести оценку параметров уравнения регрессии по формуле (9.6) с указанием нижнего индекса, относящегося к данной функции. Например, если имеется функций отклика, то оценка вектора параметров по методу наименьших квадратов может быть найдена как (9.9) Общая оценка вектора параметров многомерной функции отклика будет определяться в виде (9.10) где — символ транспонирования. Практическая часть 1. Оценка параметров одномерной функции отклика При создании моделирующего алгоритма оценивания параметров одномерной функции отклика принимаем во внимание, что наблюдения функции отклика представляют собой случайные величины, распределенные по нормальному закону [ 1 , 12 ]. При этом ранг матрицы планирования будем считать равным числу оцениваемых параметров. Входные величины в принципе могут быть произвольными. Моделирующий алгоритм будет следующим: 1. Задание базы выхода. 2. Задание входных воздействий по произвольному вероятностному закону. 3. Определение вектора выхода для r опытов. 4. Определение матрицы входа для r опытов. 5. Определение оценки вектора параметров. 6. Проверка заданной погрешности. 7. Расчет остаточной суммы квадратов. Пример 1. По заданной базе выхода (например, 6.78) и трем входным воздействий, заданных по вероятностному закону, определите параметры линейной регрессионной зависимости. Произведите оценку погрешности и рассчитайте остаточную сумму квадратов. Программный код решения примера: clear,clc % ПАРАМЕТРЫ ДИАЛОГОВОГО ОКНА inputdlg options.Resize = 'on'; options.WindowStyle = 'normal'; options.Interpreter = 'tex'; %-------------------------------------- X = inputdlg({'\bf Ввод значения функции отклика...................'}, 'База функции отклика', 1, {'6.78'}, options); Xb = str2num(char(X)); if prod(size(Xb)) = 1 | isempty(Xb) == 1 | ... isreal(Xb) | isinf(Xb) | isnan(Xb) errordlg('Отклик функции должен быть одним действительным числом',... 'Ошибка ввода функции отклика'); return end ms = inputdlg({'\bfМатематическое ожидание нормального закона:','\bfСтандартное отклонение нормального закона: .......'},... 'Параметры',1,{'0','1'},options); m = str2num(char(ms(1))); s = str2num(char(ms(2))); if isempty(m) | prod(size(m)) = 1 | isreal(m) | isinf(m) | isnan(m) errordlg('Математическое ожидание должно быть одним действительным числом',... 'Ошибка ввода математического ожидания'); return end if isempty(s) | prod(size(s)) = 1 | isreal(s) | ... isinf(s) == 1 | s <= 0 | isnan(s) errordlg('Стандартное отклонение должно быть действительным числом больше нуля',... 'Ошибка ввода стандартного отклонения'); return end % ВВОД РАЗМЕРНОСТИ ВХОДА И ЧИСЛА ОПЫТОВ u = inputdlg({'\bfЧисло входов:...',... '\bfЧисло опытов в эксперименте:..................................'},... 'Параметры эксперимента:',1,{'3','10'},options); p = str2num(char(u)); if isempty(p) | prod(size(p(1))) = 1 | isreal(p(1)) | ... isinf(p(1)) | p(1) = round(p(1)) | p(1) <= 0 | isnan(p(1)) errordlg('Число входов должно быть натуральным числом',... 'Ошибка ввода числа входов'); return end if isempty(p(2)) | prod(size(p(2))) = 1 | isnan(p(2)) | ... isreal(p(2)) | isinf(p(2)) | p(2) = round(p(2)) errordlg('Число опытов должно быть натуральным числом',... 'Ошибка ввода числа опытов'); return end if p(2) <= p(1) + 1 errordlg('Число опытов должно быть больше числа входов на 2 единицы ',... 'Ошибка ввода числа опытов'); return end %----------------------------------------------------------------- % ЗАДАНИЕ ПОГРЕШНОСТИ ПО ВЫХОДУ if abs(Xb) > 0.001 zet = inputdlg({'\bfВеличина относительной погрешности в процентах:..................'},... 'Относительная погрешность по выходу',1,{'0.5'}, options); dz = str2num(char(zet)); if isempty(dz) | prod(size(dz)) = 1 | isreal(dz) | ... isinf(dz) == 1 | dz <= 0 | isnan(dz) errordlg('Погрешность должна быть действительным числом больше нуля...',... 'Ошибка ввода относительной погрешности'); return end else zet = inputdlg({'\bfВеличина абсолютной погрешности по выходу:...'},... sprintf('Абсолютная погрешность '),1,{'0.001'}, options); dz = str2num(char(zet)); if isempty(dz) | prod(size(dz)) = 1 | isreal(dz) == 0 | ... isinf(dz) == 1 | dz <= 0 | isnan(dz) errordlg('Погрешность должна быть действительным числом больше нуля',... 'Ошибка ввода абсолютной погрешности'); return end end %----------------------------------------------------------------- Ub = [1,randn(1,p(1))]; %------------------------------------------------------------------ k = p(2); % ЧИСЛО ОПЫТОВ d = 100; n = 0; % dz - ПОГРЕШНОСТЬ ПО ВЫХОДУ while d > dz n = n + 1; % ФОРМИРОВАНИЕ ВЕКТОРА ВЫХОДА С ПОВТОРНЫМИ ОПЫТАМИ for I = 1:k R(I) = normrnd(m,s); end X = Xb + R'; %------------------------------------ U = Ub(2) + randn(k,1); % p(1) - ЧИСЛО ВХОДОВ for J = 1:p(1) - 1 U = [U,Ub(J) + rand(k,1)]; end U1 = [ones(length(U(:,1)),1),U]; %------------------------------------ b = inv(U1'*U1)*U1'*X; %------------------------------------ % РАСЧЕТНОЕ ЗНАЧЕНИЕ ВЫХОДА Xpac = Ub*b; %------------------------------------ % ОТНОСИТЕЛЬНАЯ ПОГРЕШНОСТЬ В ПРОЦЕНТАХ if abs(Xb) > 0.001 d = abs(Xpac-Xb)/abs(Xb)*100; else d = abs(Xpac-Xb); end end %------------------------------------ % ОСТАТОЧНАЯ СУММА КВАДРАТОВ Q = (X - U1*b)'*(X - U1*b); %------------------------------------ disp(sprintf('\n ОБРАБОТКА РЕЗУЛЬТАТОВ ПАССИВНОГО ЭКСПЕРИМЕНТА')); disp(sprintf(' Базовый выход Xb = %g',Xb)) fprintf(' Параметры нормального закона для формирования повторных опытов выхода:\n'); fprintf('\t математическое ожидание m = %g\n',m); fprintf('\t cтандартное отклонение s = %g\n',s); disp(sprintf('\t Число входов m = %g',length(Ub) - 1)); disp(sprintf('\t Значения входного воздействия:')); for J = 2 : length(Ub) fprintf('\t %-g\n', Ub(J)); end fprintf('\n'); if abs(Xb) > 0.001 fprintf('\t Заданная относительная погрешность в процентах: %g%s\n',dz,'%') else fprintf('\t Заданная абсолютная погрешность: %g\n',dz); end disp('--------------------------------------------------------------') fprintf('\t Число повторных опытов в эксперименте: %d\n',p(2)) fprintf('\n\t Расчетная матрица входных воздействий с повторными опытами:\n'); disp(U1) fprintf('\t Расчетный вектор выхода с повторными опытами:\n') fprintf('\t %-g\n',X) if abs(Xb) > 0.001 disp('--------------------------------------------------------------') disp(sprintf('\t Расчетный выход Xpac = %g\n\t Погрешность расчета в процентах d = %g%s',... Xpac,d,'%')) elseif abs(Xb) <= 0.001 disp('--------------------------------------------------------------') disp(sprintf('\t Расчетный выход Xpac = %g\n\t Абсолютная погрешность расчета d = %g', Xpac,d)) end fprintf('\t РАСЧЕТНЫЕ КОЭФФИЦИЕНТЫ УРАВНЕНИЯ РЕГРЕССИИ:'); for J = 1 : length(b) fprintf('\n\t b%d = %-g', J-1, b(J)); end fprintf('\n'); disp(sprintf('\t Остаточная сумма квадратов Q = %g',Q)) disp(sprintf('\t Число итераций по условию проверки погрешности n = %g',n)) fprintf('\n\t 1) Символьное уравнение регрессии:\n'); bn = length(b); fprintf('\t\t y = '),fprintf('b0 + '); for J = 1:bn - 2 fprintf('b%d*x%d + ',J,J); end fprintf('b%d*x%d\n',bn - 1,bn - 1); fprintf('\t 2) Уравнение регрессии с числовыми оценочными параметрами:\n') bn = length(b); fprintf(' y = '),fprintf('%g + ',b(1)); for J = 2 : (bn - 1) fprintf('(%g)*x%d + ',b(J),J - 1); end fprintf('(%g)*x%d\n',b(end), bn - 1); disp('--------------------------------------------------------'); Ниже приводится одна из возможных реализаций выполнения программы (т. к. выход и входы смоделированы по случайным законам). ОБРАБОТКА РЕЗУЛЬТАТОВ ПАССИВНОГО ЭКСПЕРИМЕНТА Базовый выход: Xb = 6.78 Параметры нормального закона для формирования повторных опытов выхода: математическое ожидание: m = 0 cтандартное отклонение: s = 1 Число входов: m = 3 Значения входного воздействия: -1.01743 -1.4031 -0.260899 Заданная относительная погрешность в процентах: 0.5% -------------------------------------------------------------- Число повторных опытов в эксперименте: 10 Расчетная матрица входных воздействий с повторными опытами: 1.0000 -0.0259 1.8706 -0.8428 1.0000 -2.0962 1.5201 -0.2298 1.0000 -0.7237 1.4284 -0.1887 1.0000 -1.4151 1.6469 -0.7250 1.0000 -0.8569 1.3781 -0.1382 1.0000 -0.3895 1.9245 -0.0297 1.0000 -2.4837 1.0058 -0.9198 1.0000 0.0256 1.7664 -0.6273 1.0000 -1.1517 1.6920 -0.8416 1.0000 -0.9158 1.1305 -0.5560 Расчетный вектор выхода с повторными опытами: 5.86991 5.30057 9.27119 9.05969 6.35003 9.75507 8.46856 6.92237 9.01667 5.847 -------------------------------------------------------------- Расчетный выход: Xpac = 6.80021 Погрешность расчета в процентах: d = 0.298114% РАСЧЕТНЫЕ КОЭФФИЦИЕНТЫ УРАВНЕНИЯ РЕГРЕССИИ: b0 = 6.48866 b1 = -0.159952 b2 = -0.0359629 b3 = -0.376992 Остаточная сумма квадратов: Q = 9.89971 Число итераций по условию проверки погрешности: n = 184 1) Символьное уравнение регрессии: y = b0 + b1*x1 + b2*x2 + b3*x3 2) Уравнение регрессии с числовыми оценочными параметрами: y = 6.48866 + (-0.159952)*x1 + (-0.0359629)*x2 + (-0.376992)*x3 -------------------------------------------------------------- Задание 1 1. В программе объясните цикл проверки условия погрешности по выходу. 2. Проверьте программу для случая задания абсолютной погрешности по выходу. 3. Видоизмените программу так, чтобы расчеты производились без формирования столбца единиц в матрице входных воздействий с повторными опытами. 4. Расчет вектора оценочных параметров произведите по встроенной функции системы MATLAB (см. ). 5. Расчет вектора оценочных параметров произведите по встроенной функции системы MATLAB (см. ). Для функции нет необходимости формировать столбец фиктивной переменной, равной единице. 6. Расчет вектора оценочных параметров произведите по встроенной функции системы MATLAB (см. ). Произвести также расчет значений функции отклика с помощью функции системы MATLAB (см. ). Примечание. Описание функций , , на русском языке можно найти на сайтессылка: http://matlab.exponenta.ru - http://matlab.exponenta.ruв разделе Statistics Toolbox. 2. Аппроксимация функций степенными полиномами С помощью метода наименьших квадратов возможно определять параметры уравнения кривых, заданных некоторыми своими точками, при условии, что вид уравнения известен, а неизвестными выступают параметры этого уравнения. При этом всякая аналитическая функция может быть приближена полиномом. Поэтому для подгонки нелинейных данных возможно использовать построение полинома методом наименьших квадратов. Для более точной подгонки степень полинома приходится увеличивать. Правда, если данные не проявляют полиномиальной природы, то в результате получается кривая, которая будет сильно осциллировать [15] . Этот феномен называется полиномиальное раскачивание [15 , с. 301]. Оно явно наблюдается у полиномов высокой степени. Из этих соображений полиномы степени 6 или выше редко используются, если известно, что истинная функция, выражающая зависимость, является полиномом. Но в то же время приближение данных полиномами служит основным инструментом для случая, когда уравнение кривой неизвестно совсем. Тогда в задачу входит определение коэффициентов (параметров) выбранного полинома, которым осуществляется подгонка. В любом случае метод наименьших квадратов дает достаточно приемлемые результаты. Более того, когда данные представляют собой некоторую "россыпь", с помощью метода наименьших квадратов возможно провести кривую между этими данными, которые наилучшим образом минимизируют квадрат невязки по этим данным. Пример 2. Для функции вида известны несколько значений точек на плоскости. По этим значениям рассчитайте параметры (коэффициенты) заданной функции. Программный код решения примера: clear,clc try global hep delete(hep); end options.Resize = 'on'; options.WindowStyle = 'normal'; options.Interpreter = 'tex'; f = inputdlg({sprintf('\\bfВектор базисных функций:'),... sprintf('\\bfОбласть определения:'),... sprintf('\\bfЧисловые значения параметров b0,b1,...'),... sprintf('\\bfЗначения аргумента (не менее 4 значений):...')},... 'ИСХОДНЫЕ ДАННЫЕ',1,... {'x;x^3;sqrt(x)','[0 2]','-5;1.2;3.2;-12.4','0,0.12,0.7,1.34,0.45,1.678 '},options) if isempty(f) errordlg('Возможно, нажата клавиша ''Cansel'' ',... 'Ошибка'); return end y1 = sym(char(f(1))); x = str2num(char(f(2))); b = str2num(char(f(3))); es = str2num(char(f(4))); Y2 = char(y1); q = find(Y2 == 'x'); if length(q) + 1 = length(str2num(char(f(3))) ) errordlg('Число введенных параметров не соответствует числу параметров введенной функции',... 'Ошибка ввода параметров'); return end %---------------- Общее число точек --------------------- if length(es) == 1 if es < length(b) errordlg(sprintf('Общее число точек должно быть не менее числа %d - числа параметров',length(b)),... 'Ошибка ввода числа точек функции'); return end fprintf('\n ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ЗАДАННОЙ ФУНКЦИИ\n'); fprintf('\t Вид параметрических функций:\n\t %s\n',Y2); fprintf('\t Область определения функции: [%g %g]\n',x(1),x(2)); fprintf('\t Значения введенных параметров:\n'); for J = 1:length(b) fprintf('\t%g\t',b(J)'); end fprintf('\n'); disp('-----------------------------------------------------') er1 = inline(sprintf('[1,%s]',char(f(1)))); erV1 = vectorize(inline(sprintf('[1,%s]',char(f(1))))); erV = vectorize(inline(sprintf('[%s]',char(f(1))))); con = char(erV); n = find(con == ';'); con(n) = ','; con1 = char(er1); n1 = find(con1 == ';'); con1(n1) = ','; symV = sym(con1); wone = inline(con); V =x(1) + (x(2) - x(1))*rand(es,1); fprintf('\t Значения аргумента и параметрических функций:\n'); X = [ones(length(V),1),wone(V)]; disp([V,X]) XX = []; for I = 1:length(V) XX = [XX;b'.*X(I,:)]; end Y = sum(XX,2); fprintf('\t Значения аргумента и заданной функции:\n'); disp([V,Y]) B = regress(Y,X); fprintf('\t Расчетные значения параметров заданной функции:\n'); for J = 1:length(B) fprintf('\t b%d = %g\n',J - 1, B(J)); end y = []; for J = 1:length(symV) if J < length(symV) if J == 1 y = [y,sprintf('b%d%s + ',J - 1,'')]; else y = [y,sprintf('b%d*%s + ',J - 1,char(symV(J)))]; end elseif J == length(symV) y = [y,sprintf('b%d*%s',J - 1,char(symV(end)))]; end end fprintf('\n\t Заданный вид функциональной зависимости:\n\t y = %s\n',y) helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'),' ') %----------------- Несколько точек ------------------------------ elseif length(es) >= length(b) if min(es) < x(1) | max(es) > x(2) errordlg(sprintf('Значения аргумента функции выходят за область определения [%g, %g]',x(1),x(2)),... 'Ошибка ввода значений аргумента функции'); return end fprintf('\n ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ЗАДАННОЙ ФУНКЦИИ\n'); fprintf(' Вид параметрических функций:\n\t %s\n',Y2); fprintf(' Область определения функции: [%g %g]\n',x(1),x(2)); disp('-------------------------------------------------------------') er1 = inline(sprintf('[1,%s]',char(f(1)))); erV1 = vectorize(inline(sprintf('[1,%s]',char(f(1))))); erV = vectorize(inline(sprintf('[%s]',char(f(1))))); con = char(erV); n = find(con == ';'); con(n) = ','; con1 = char(er1); n1 = find(con1 == ';'); con1(n1) = ','; symV = sym(con1); wone = inline(con); V = es'; fprintf(' Значения аргумента и значения параметрических функций:\n'); X = [ones(length(V),1),wone(V)]; disp([V,X]) XX = []; for I = 1:length(V) XX = [XX;b'.*X(I,:)]; end Y = sum(XX,2); fprintf(' Значения аргумента и значения заданной функции:\n'); disp([V,Y]) B = regress(Y,X); fprintf(' Расчетные значения параметров заданной функции:\n'); for J = 1:length(B) fprintf('\t b%d = %g\n',J - 1,B(J)); end y = []; for J = 1:length(symV) if J < length(symV) if J == 1 y = [y,sprintf('b%d%s + ',J - 1,'')]; else y = [y,sprintf('b%d*%s + ',J - 1,char(symV(J)))]; end elseif J == length(symV) y = [y,sprintf('b%d*%s',J - 1,char(symV(end)))]; end end fprintf('\n Заданный вид функциональной зависимости:\n\t y = %s\n',y); global hep hep = helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'), ''); else errordlg(sprintf('Необходимо общее число аргументов не менее %d',length(b)), 'Ошибка!') return end Задание 2 1. Измените область определения параметрических функций, например, 2. Постройте диаграмму заданной функциональной зависимости. 3. В программе задайте "свои" параметрические функции. 4. Определение параметров (коэффициентов) заданной функции произведите через решение нормального уравнения с помощью метода исключения Гаусса. Пример: 5. В программе предусмотрите "зашумление" значений заданной функции по нормальному и равномерному вероятностным законам. Рассчитайте значения параметров (коэффициентов) заданной функции и подсчитайте относительную погрешность расчета коэффициентов. Подсчитайте также остаточную сумму квадратов; 6. Изучите функции , по справочной системе пакета MATLAB. Пример 3. Функциональная экспериментальная зависимость задана массивом точек (значений). Произвести аппроксимацию данной функциональной зависимости степенным полиномом. Постройте графики функций. Значения экспериментальной зависимости: x f(x) 0; 149.0; 5; 146.3; 10; 145.8; 15; 149.8; 20; 154.9; 25; 164.7; 30; 178.2; 35; 196.2; 40; 221.8; 45; 259.8; 50; 299.2; 55; 331.5; 60; 338.1; 65; 304.4; 70; 238.5; 72; 205.3; 74; 170.9; 75; 153.9; 76; 136.6; 78; 103.9; 80; 74.2; 82; 48.6; 84; 29.7; 85; 19.9; 86; 12.4; 88; 3.1; 90; 0; Решение примера выполним с помощью стандартных функций MATLAB, таких как , , , [8] Программный код решения примера: clear,clc,close all x = ... [0;5;10;15;20;25;30;35;40;45;50;55;60;65;70;72;74;75;76;78;80;82;84;85;86;88;90]; f = ... [149.0;146.3;145.8;149.8;154.9;164.7;178.2;196.2;221.8;259.8;299.2;331.5;338.1;304.4;238.5;205.3;170.9;153.9; 136.6;103.9;74.2;48.6;29.7;19.9;12.4;3.1;0]; n = length(x); %% length(x) == length(f) Lmin = f'*f; %% Квадратичный критерий - сумма квадратов line(x,f,'linew', 2) grid on xlabel('\bf - - - - - - - - - x - - - - - - - - -'); ylabel('\bf f(x)'); qprob = 0.4; % Уровень значимости %%%% Выбор степени полинома for m = 1 : n psi(:, m) = x.^(m-1); % psi - Специальная функция (полигамма) if m > 1 for k = 1 : m-1 psi(:, m) = psi(:, m) - ((psi(:,m))'*psi(:,k))*psi(:, k); end end psi(:, m) = psi(:,m)/norm(psi(:,m)); % нормировка bk = f'*psi(:,m); DOLD = Lmin/(n-m); %%%Старая дисперсия Lmin = Lmin - bk^2; Dnew = Lmin/(n-m-1); %%% Новая дисперсия % fprintf('\n'); if m > 1 % проверка степени начиная с 2 if Dnew/DOLD < finv(qprob, n-m, n-m-1) % fprintf(' Степень %d следует учитывать\n', m-1); stp = m-1; else % fprintf(' Степень %d не нужно учитывать\n', m-1); if m > 2 stp = m-2; break; else stp = m - 1; break; end end end end PO = polyfit(x, f, stp); fprintf('\n Коэффициенты аппроксимируещего полинома:\n'); for J = 1 : length(PO) fprintf('%15g\n', PO(J)); end fprintf('\n'); yPO = polyval(PO, x); line(x, yPO, 'marker', 'o', 'color','r') title(sprintf('%s Аппроксимация экспериментальной кривой полиномом %d степени','\bf', stp)); legend('\bfЭксперимент','\bfАппроксимация', 'location', 'best') Диаграммы заданной экспериментальной кривой и аппроксимирующего полинома показаны на рис. 9.2 Рис. 9.2. Результат аппроксимации кривой полиномом Задание 3 1. Проверьте результат выполнения программы при уменьшении и увеличении уровня значимости. 2. Постройте аппроксимирующую функцию непосредственно по аппроксимирующему полиному (без функции ). 3. Запишите в текстовый файл полученный аппроксимирующий полином. Имя файла задайте в виде compX.txt, где — номер компьютера, за которым выполняется лабораторная работа (1, 2, 3, ...). 4. Создайте массив значений функции на отрезке и произведите аппроксимацию полиномами следующих порядков: 2-го, 3-го, 4-го. Рассчитайте и постройте приведенную погрешность аппроксимации. Примечание. В качестве приведенной погрешности примите разность в абсо-лютных значениях между аппроксимирующей функцией и истинной функ-ции, отнесенной к максимуму модуля истинной функции (например, для функции ). 3. Оценка параметров многомерной функции ОТ-клика Пример 4. Поставим задачу оценки параметров , , , многомерной функции отклика вида по известным средним значениям независимых аргументов , число которых , и по известным средним значениям самих функций , число которых . Будем считать, что совокупность повторных значений функций распределена по нормальному закону, а аргументы функций могут быть распределены произвольно. Параметры многомерной функции отклика образуют матрицу размера . Эту матрицу обозначим как . Как видно, функции отклика являются линейными как по своим аргумента, так и по параметрам. В связи с этим для оценки параметров по измерениям аргументов и функций отклика (повторным значениям) целесообразно применить метод наименьших квадратов. Программный код решения примера: clear,clc %--------------------------------------------------------------- options.Resize = 'on'; options.WindowStyle = 'normal'; options.Interpreter = 'tex'; yx = inputdlg({'\bfСредние значения функций отклика:',... '\bfСредние значения аргументов функции отклика:',... '\bfЧисло измерений:',... '\bfВеличина стандартного отклонения функции откли-ка:.....'},... 'ИСХОДНЫЕ ДАННЫЕ',1,{'1,2,3,4,5,6,7,8','11,23,4,7,9','50','0.1'},options); % БАЗОВЫЕ ВЕКТОРЫ ВЫХОДА, ВХОДА, ЧИСЛО ИЗМЕРЕНИЙ, % СТАНДАРТНОЕ ОТКЛОНЕНИЕ Yb = str2num(char(yx(1))); Ub = str2num(char(yx(2))); r = str2num(char(yx(3))); s = str2num(char(yx(4))); %------------------------------------------------------------------- % ПРОВЕРКА ВВОДА ДАННЫХ if ( isempty(Yb) == 1 | isreal(Yb) == 0 | isinf(Yb) == 1 ) errordlg(sprintf('Значения функций должны составлять вектор действительных чисел'),'Ошибка ввода значений функций') return elseif ( min(size(Yb)) > 1 & max(size(Yb)) > 1 ) errordlg(sprintf('Значения функций должны составлять вектор, не матрицу'),'Ошибка ввода значений функций') return end %---------------------------------------------------------------- %% ПРИВЕДЕНИЕ К ВЕКТОРУ-СТРОКЕ [ry,vy] = size(Yb); if ( ry == 1 ) Yb = Yb; elseif ( ry > 1 ) Yb = Yb'; end if ( isempty(Ub) | isreal(Ub) | isinf(Ub) ) errordlg(sprintf('Значения аргументов должны составлять вектор действительных чисел'), ... 'Ошибка ввода значений аргументов') return elseif ( min(size(Ub)) > 1 & max(size(Ub)) > 1 ) errordlg('Значения аргументов функций должны составлять вектор, не матрицу','Ошибка ввода значений аргументов функций') return elseif ( isempty(find(Ub == 0)) == 0 ) errordlg(sprintf('Значения аргументов функций не должны быть нулевыми'),'Ошибка ввода значений аргументов') return end %-------------------------------------------------------------- %% ПРИВЕДЕНИЕ К ВЕКТОРУ-СТРОКЕ [ru,vu] = size(Ub); if ( ru == 1 ) Ub = Ub; elseif ( ru > 1 ) Ub = Ub'; end if ( isempty(r) | isreal(r) | isinf(r) | r <= 0 ) errordlg('Число измерений должно быть натуральным числом',... 'Ошибка ввода числа измерений') return elseif ( length(r) > 1 | r = round(r)) errordlg('Число измерений должно быть одним натуральным числом',... 'Ошибка ввода числа измерений') return elseif ( r <= length(Yb)*(length(Ub) + 1) ) errordlg(sprintf('Число измерений должно быть не меньше произведения числа функций на число аргументов плюс 1: больше %d',... length(Yb)*(length(Ub) + 1)), 'Ошибка ввода числа измерений') return end %----------------------------------------------------------------- if ( isempty(s) | isreal(s) | isinf(s) ) errordlg('Стандартное отклонение должно быть действительным числом',... 'Ошибка ввода стандартного отклонения') return elseif ( length(s) > 1 ) errordlg('Стандартное отклонение должно быть одним действительным числом','Ошибка ввода стандартного отклонения') return elseif ( s <= 0 ) errordlg('Стандартное отклонение должно быть больше нуля',... 'Ошибка ввода отклонения') return end %---------------------------------------------------------------- % % ЗАДАНИЕ r ИЗМЕРЕНИЙ ФУНКЦИЙ ОТКЛИКА Y = []; for k = 1 : length(Yb) for I = 1 : r Y(:,k) = Yb(k) + normrnd(0,s,r,1); end end % % ЗАДАНИЕ r ИЗМЕРЕНИЙ АРГУМЕНТОВ ФУНКЦИИ ОТКЛИКА X = []; for m = 1 : length(Ub) for J = 1 : r X(:,m) = Ub(m) + rand(r,1); end end %% ВВЕДЕНИЕ ФИКТИВНОГО АРГУМЕНТА, РАВНОГО ЕДИНИЦЕ X1 = [ones(r,1),X]; %% РАСЧЕТ ПАРАМЕТРОВ ФУНКЦИИ ОТКЛИКА B = (inv(X1'*X1)*X1'*Y)'; %% ВЫВОД РЕЗУЛЬТАТОВ fprintf('\n\t ИСХОДНЫЕ ДАННЫЕ ЭКСПЕРИМЕНТА:\n') fprintf('\t Вектор-строка значений функции отклика:\n') fprintf('\t %g',Yb) fprintf('\n') fprintf('\t Вектор-строка значений аргументов:\n') fprintf('\t %g',Ub) fprintf('\n') fprintf('\t Заданное число повторных измерений: r = %d\n',r) fprintf('\t Стандартное отклонение повторных измерений функций: s = %g\n',s) disp('-----------------------------------------------------------') fprintf('\t РЕЗУЛЬТАТЫ ПАССИВНОГО МНОГОМЕРНОГО ЭКСПЕРИМЕНТА\n') fprintf('\t ОЦЕНКА МАТРИЦЫ ПАРАМЕТРОВ:\n') disp(B) %% РАСЧЕТНЫЕ ЗНАЧЕНИЯ ФУНКЦИЙ Ybpac = B*[1,Ub]'; disp('------------------------------------------------------------') fprintf('\t СРАВНЕНИЕ ЗНАЧЕНИЙ ЗАДАННЫХ ФУНКЦИЙ И РАСЧЕТНЫХ:\n') disp([Yb' Ybpac]) disp('------------------------------------------------------------') helpdlg(sprintf('Результаты выполнения программы смотрите в командном окне'),'') Пример выполнения программы без демонстрации диалоговых окон ИСХОДНЫЕ ДАННЫЕ ЭКСПЕРИМЕНТА Вектор-строка значений функции отклика: 1 2 3 4 5 6 7 8 Вектор-строка значений аргументов: 11 23 4 7 9 Заданное число повторных измерений: r = 50 Стандартное отклонение повторных измерений функций: s = 0.1 ----------------------------------------------------------- РЕЗУЛЬТАТЫ ПАССИВНОГО МНОГОМЕРНОГО ЭКСПЕРИМЕНТА ОЦЕНКА МАТРИЦЫ ПАРАМЕТРОВ: -1.5280 -0.0296 0.0969 -0.0103 -0.0108 0.0767 3.1418 -0.0457 -0.0015 0.0794 -0.0597 -0.0519 2.2212 -0.0010 0.0370 -0.0022 0.0191 -0.0225 4.3340 -0.1121 0.0989 -0.0480 0.0087 -0.1287 1.3881 0.0194 0.0830 -0.0053 0.0519 0.1112 8.1169 -0.0406 -0.1085 0.0382 -0.0001 0.0760 9.9112 -0.0201 -0.0095 0.0651 -0.0843 -0.0143 8.1976 0.0707 -0.0004 -0.0371 -0.0408 -0.0563 ------------------------------------------------------------ СРАВНЕНИЕ ЗНАЧЕНИЙ ЗАДАННЫХ ФУНКЦИЙ И РАСЧЕТНЫХ: 1.0000 0.9483 2.0000 2.0366 3.0000 2.9832 4.0000 4.0865 5.0000 4.8532 6.0000 6.0105 9.0000 9.0124 8.0000 8.0252 ------------------------------------------------------------ Примечание. Обратите внимание на зависимость полученных результатов относительно задаваемого стандартного отклонения, применяемого для формирования повторных измерений функций отклика. Задание 4 1. Результат выполнения программы запишите в текстовый файл с именем compX.txt, где — номер компьютера, за которым выполняется лабораторная работа (1, 2, 3, ...). 2. Дополните программу вычислением отдельных составляющих заданной многомерной функции отклика. 3. Напишите фрагмент программы по расчету относительной и, возможно, абсолютной погрешности относительно заданных значений функций отклика. 4. Рассчитайте остаточную сумму квадратов по составляющим функции отклика. 5. Видоизмените программу для оценки параметров многомерной функции отклика для случая, когда в уравнениях отсутствуют свободные слагаемые Контрольные вопросы 1. В чем заключается основная задача пассивного эксперимента? |