1. Что такое функция отклика?
2. В чем различие между пассивным и активным экспериментами?
3. Какие задачи решает активный эксперимент?
4. Что называется градиентом функции отклика?
5. В чем заключается методика поиска экстремума функции отклика?
6. Почему используются линейные модели при поиске экстремума
функции отклика?
7. Как осуществляется построение факторного эксперимента в окрестности заданной точки факторного пространства?
8. На каком этапе поиска экстремума функции отклика применяется метод наименьших квадратов?
Оценивание параметров линейной модели по наблюдениям неполного ранга
Цель работы: изучить способы оценки параметров линейной регрессионной модели в случае вырожденной информационной матрицы нормального уравнения. Среда программирования —
MATLAB.
Теоретическая часть
Пусть линейная модель наблюдений имеет следующий вид в векторно-
матричной форме:
(11.1)
где:
— -мерный вектор наблюдений;
— матрица известных коэффициентов порядка
(матрица регрессоров);
— -мерный вектор неизвестных параметров, подлежащих определению;
— -мерный случайный вектор (помехи);
— неизвестный параметр;
— единичная матрица;
— операторы взятия математического ожидания и дисперсии
[1]
В случае применения метода наименьших квадратов для оценки параметров модели (11.1) приходим к нормальному уравнению
(11.2)
Если
, то информационная матрица является вырожденной (т. е. ее детерминант равен нулю) и множество решений нормального уравнения (11.2) будет бесконечным
[1]
. В этом случае модель наблюдений (11.1) называют моделью наблюдений неполного
ранга.
Оценка общего решения нормального уравнения с вырожденной информационной матрицей может быть получена с помощью обобщенной обратной матрицы, которую обозначают как и
называют -обратной матрицей. Обобщенная обратная матрица всегда существует и не обязательно единственная [
1
,
17
].
Определение. Пусть
-матрица произвольного ранга.
Обобщенной обратной (или -обратной) матрицей для называется матрица порядка
, такая, что
(11.3)
Нормальное уравнение (11.2) запишем в виде
(11.4)
где:
— квадратная матрица порядка ;
— -мерный вектор.
Если
, то матрица — вырожденная.
Общее решение неоднородного уравнения (11.4) есть сумма любого
частного решения этого уравнения и общего решения однородного
уравнения
(11.5)
Известно, что если — обобщенная матрица для матрицы , то частное решение уравнения (11.4) имеет вид (11.6) Общее решение однородного уравнения (11.5) имеет вид (11.7) где: — идемпотентная матрица ( ); — единичная матрица порядка ; — произвольный вещественный вектор порядка , в том числе и равный тождественно нулю [1] Общее решение уравнения (11.4) есть сумма частного решения (11.6) и решения однородного уравнения (11.7), т. е. (11.8) Алгоритм построения общего решения нормального уравнения и обобщенной обратной матрицы S 1. В матрице неполного ранга определяют линейно независимые столбцы. Пусть это будут первые -столбцов. 2. Матрицу представляют в виде блочной матрицы с подматрицами и , т. е. , где — матрица размера 3. Матрицу представляют в виде 4. Обобщенную обратную матрицу определяют в виде
5. Общее решение нормального уравнения — оценку вектора
параметров определяют по формуле (11.8) или в следующем виде:
где
— единичная матрица порядка
6. Вычисляют остаточную сумму квадратов RSS (Residual Sum of
Squares):
Практическая часть
2.1. Оценка параметров модели наблюдений неполного ранга с помощью обобщенной обратной матрицы
Пример 1. Для заданной матрицы регрессоров неполного ранга и вектора значений функции отклика произведите оценку линейной регрессионной модели. Рассчитайте остаточную сумму квадратов RSS.
Постройте зависимость RSS от одного из компонент вектора коэффициентов (параметров) линейной регрессионной модели.
Исходные данные
Матрица регрессоров:
вектор значений функции отклика (вектор наблюдений):
Сначала следует определить линейно зависимые столбцы и произвести преобразование матрицы с целью расположения первых линейно независимых столбцов. Величина будет равняться рангу матрицы.
Как известно, столбцы матрицы — линейно зависимые, если значения одного столбца отличаются от значений другого столбца на постоянный множитель. Порядок расположения в матрице регрессоров не отразится на результате произведения этой матрицы на искомые коэффициенты
линейной модели, если в таком же порядке произвести перестановку искомых коэффициентов.
Программный код решения примера:
clear,clc,close
% Исходные данные
% Матрица регрессоров
X = [1.340, 1.630, -4.240, 2.412, 5.088;
1.610, 1.900, -3.970, 2.898, 4.764;
1.470, 1.760, -4.110, 2.646, 4.932;
2.610, 2.900, -2.970, 4.698, 3.564;
3.000, 3.290, -2.580, 5.400, 3.096;
2.040, 2.330, -3.540, 3.672, 4.248;
1.060, 1.470, -4.520, 1.908, 5.424];
% Вектор наблюдений - вектор выхода
Y = [13.784;13.617;12.968;14.713;13.864;12.847;14.420];
% Размерность вектора наблюдений n = length(Y);
% Размерность искомого вектора коэффициентов линейной модели p = size(X,2);
%%-------------------------------------------
% Вычисление ранга матрицы регрессоров rankX = rank(X);
if rankX == p %% модель полного ранга fprintf('\n Размерность вектора параметров модели: p = %d\n Модель наблюдений полного ранга: rank(X) = %d\n', p, rankX);
B = inv(X'*X)*X'*Y;
disp(' ')
disp(' Оценка вектора параметров B: ')
disp(B)
YY = X*B;
fprintf('\n Y - заданный вектор наблюдений;\n');
fprintf(' YY - расчетный вектор наблюдений;\n');
fprintf(' Error - относительная погрешность;\n');
fprintf('\tY\t\t\tYY\t\t Error\n');
erd = abs(Y-YY)./abs(Y)*100;
for J = 1 : length(Y)
fprintf(' %1.4f\t %1.4f\t %1.4f%%\n', Y(J), YY(J), erd(J));
end
RESS = (Y-YY)'*(Y-YY);
Cp = mean(abs(Y-YY)./abs(Y)*100);
fprintf('\n Величина RSS: %g\n', RESS);
fprintf(' Усредненная погрешность расчета вектора наблюдения: %g%%\n', Cp);
return
end
%% модель наблюдений неполного ранга fprintf('\n Размерность вектора параметров модели: p = %d\n',p);
fprintf(' Модель наблюдений неполного ранга: rank(X) = %d\n', rankX);
% Определение номеров линейно зависимых столбцов k1 = 0;
for N = 1 : p for J = 1 : p op = abs(X (:,J))./abs(X (:,N));
Pr(1:n,1) = op;
str = sprintf('%1.5f', Pr(1));
str2 = sprintf('%1.5f', Pr(end));
if str2num(str) == str2num(str2) & str2num(str) = 1 k1 = k1 + 1;
NJ(k1,1:2) = [N,J];
end end end k2 = 0;
en = size(NJ,1);
for I = 1 : size(NJ,1)
for J = 1 : size(NJ,2)
if NJ(I,J) == NJ(en-J+1,2)
k2 = k2 + 1;
EX(k2) = (en-J+1);
end end end
NJ(EX,:) = []
for J = 1 : size(NJ,1)
fprintf(' Столбец %d линейно зависит от столбца %d\n', NJ(J,2), NJ(J,1));
end
%%-------------------------------------------
X0 = X(:,1:rank(X));
X1 = X(:, rank(X)+1:end);
S = X'*X; %% информационная матрица
% Формирование g-обратной матрицы
B11 = inv(X0'*X0);
B12 = zeros(rank(X), p-rank(X));
B21 = zeros(p-rank(X), rank(X));
B22 = zeros(p-rank(X),p -rank(X) );
Sg = [B11, B12;B21,B22] % g-обратная матрица
H = Sg*S;
z = randn(size(X,2), 1); %%случайный вещественный вектор
P11 = inv(X0'*X0)*X0'*Y;
P21 = zeros(p - rank(X), 1);
% Z11 = zeros(rank(X));
% Z12 = inv(X0'*X0)*X0'*X1;
% Z21 = zeros(p - rank(X), rank(X));
% Z22 = -eye((p - rank(X)));
% B = [P11;P21]+ [Z11,Z12;Z21,Z22]*z
Ep = eye(p);
B = [P11;P21] + [H - Ep]*z;
disp(' ')
disp(' Оценка вектора параметров B: ')
disp(B)
YY = X*B;
fprintf('\n Y - заданный вектор наблюдений;\n');
fprintf(' YY - расчетный вектор наблюдений;\n');
fprintf(' Error - относительная погрешность;\n');
fprintf('\tY\t\t\tYY\t\t Error\n');
erd = abs(Y-YY)./abs(Y)*100;
for J = 1 : length(Y)
fprintf(' %1.4f\t %1.4f\t %1.4f%%\n', Y(J), YY(J), erd(J));
end
RESS = (Y-YY)'*(Y-YY);
Cp = mean(abs(Y-YY)./abs(Y)*100);
fprintf('\n Величина RSS: %g\n', RESS);
fprintf(' Усредненная погрешность расчета вектора наблюдения: %g%%\n', Cp);
Bopt = B;
k = 2;
d = 0;
for J = 0.79*Bopt(k) : 0.01 : 1.21*Bopt(k)
d = d + 1;
B(k) = J;
YY = X*B;
RESS(d) = (Y-YY)'*(Y-YY);
Xres(d) = J;
end kmin = find(RESS == min(RESS));
bmin = Xres(kmin);
fig1 = figure(1);
set(fig1,'name','Зависимость RSS от параметра')
line(Xres, RESS, 'linew', 2)
grid on ch = sprintf('%sОптимальное значение %d-го параметра: %s(%d) = %g','\bf\fontsize{11}\fontname{times}', k, strcat('\beta', '_{opt}'), k,bmin);
title(ch);
ylabel('\bf\fontsize{12}\fontname{times}RSS');
ch2=sprintf('%s_%d','\bf\fontsize{12}\fontname{times}\beta',k);
xlabel(ch2);
В программе предусмотрен вариант расчета вектора параметров линейной регрессионной модели для случая наблюдений полного ранга.
Если имеется полный ранг, то расчет ведется по решению нормального уравнения, в котором информационная матрица — не вырожденная.
Притом с помощью оператора прекращается дальнейшее выполнение программы.
В программе усредненная погрешность расчета вектора параметра – это средняя величина от столбца относительной погрешности в процентах.
В программе произведено изменение 2-го параметра регрессионной модели. Его изменение выполнено слева и справа от оптимального значения.
Результат выполнения программы
Размерность вектора параметров модели: p = 5
Модель наблюдений неполного ранга: rank(X) = 3
Столбец 4 линейно зависит от столбца 1
Столбец 5 линейно зависит от столбца 3
Sg =
116.6611 -112.7327 -7.3650 0 0
-112.7327 109.0775 7.1933 0 0
-7.3650 7.1933 0.5164 0 0 0 0 0 0 0 0 0 0 0 0
Оценка вектора параметров B:
-7.8809 10.5189
-0.9992
-0.2349 0.5978
Y - заданный вектор наблюдений;
YY - расчетный вектор наблюдений;
Error - относительная погрешность;
Y
YY
Error
13.7840 13.2973 3.5308%
13.6170 13.4319 1.3592%
12.9680 13.3621 3.0392%
14.7130 13.9305 5.3187%
13.8640 14.1249 1.8817%
12.8470 13.6463 6.2216%
14.4200 14.4200 0.0000%
Величина RSS: 1.74575
Усредненная погрешность расчета вектора наблюдения: 3.05018%
Диаграмма зависимости RSS от величины 2-го компонента вектора параметров показана на рис. 11.1
Рис. 11.1. Зависимость RSS от 2-го компонента вектора параметров модели
Задание 1
1. Вектор определите с помощью функции
(см.
). В качестве параметра функции примите номер компьютера, на котором выполняется лабораторная работа (1, 2, 3, ...).
2. Вектор определите с помощью равномерного закона распределения вероятностей из интервала
, где — номер компьютера, на котором выполняется лабораторная работа (1, 2, 3,
...).
3. При расчете вектора наблюдений (
) произведите его "зашумление" случайным вектором, распределенным по стандартному нормальному закону.
4. Постройте зависимость RSS от числа прогонов программы. Число принять от 1 до
, где — номер компьютера, на котором выполняется лабораторная работа (1, 2, 3, ...).
5. Напишите программу по определению линейно независимых столбцов и расположению их на первые мест матрицы регрессоров.
6. Найдите оценку параметров модели неполного ранга (по данным примера 1) с помощью псевдообратной матрицы Мура-Пенроуза
(см.
).
Пример 2. Вычислите коэффициенты уравнения линейной регрессии при следующем заданном множестве измерений эксперимента:
Таблица 11.1.
Результаты измерений эксперимента
0.620 0.400 0.420 0.820 0.660 0.720 0.380 0.520 0.450 0.690 0.550 12.000 14.200 14.600 12.100 10.800 8.200 13.000 10.500 8.800 17.000 14.400 5.890 3.800 3.990 7.790 6.270 6.840 3.610 4.940 4.275 6.555 5.225 51.60 49.90 48.50 50.60 49.70 48.80 42.60 45.90 37.80 64.80 53.40
Данные, приведенные в таблице 11.1,
взяты с некоторыми изменениями из
[7]
Составим матрицу регрессоров :
Составим вектор наблюдений :
Программный код решения примера:
clear,clc
%% Исходные данные
%% Матрица регрессоров
X = [0.620, 12.000, 5.890;
0.400, 14.200, 3.800;
0.420, 14.600, 3.990;
0.820, 12.100, 7.790;
0.660, 10.800, 6.270;
0.720, 8.200, 6.840;
0.380, 13.000, 3.610;
0.520, 10.500, 4.940;
0.450, 8.800, 4.275;
0.690, 17.000, 6.555;
0.550, 14.400, 5.225;
0.360, 12.800, 3.420];
%% Вектор наблюдений - вектор выхода
Y = [51.60;49.90;48.50;50.60;49.70;48.80;42.60;45.90;37.80;...
64.80;53.40;45.30];
% Размерность вектора параметров - коэффициентов модели p = size(X,2);
% Размерность вектора наблюдений n = length(Y);
rankX = rank(X);
if rankX == p %% модель полного ранга fprintf('\n Размерность вектора параметров модели: p = %d\n Модель наблюдений полного ранга: rank(X) = %d\n', p, rankX);
B = inv(X'*X)*X'*Y;
disp(' ')
disp(' Оценка вектора параметров B: ')
disp(B)
return end
%% модель наблюдений неполного ранга fprintf('\n Размерность вектора параметров модели: p = %d\n',p);
fprintf(' Модель наблюдений неполного ранга: rank(X) = %d\n', rankX);
% Определение номеров линейно зависимых столбцов k1 = 0;
for N = 1 : p for J = 1 : p op = abs(X (:,J))./abs(X (:,N));
Pr(1:n,1) = op;
str = sprintf('%1.4f', Pr(1));
str2 = sprintf('%1.4f', Pr(end));
if str2num(str) == str2num(str2) & str2num(str) = 1 k1 = k1 + 1;
NJ(k1,1:2) = [N,J];
end end end k2 = 0;
en = size(NJ,1);
for I = 1 : size(NJ,1)
for J = 1 : size(NJ,2)
if NJ(I,J) == NJ(en-J+1,2)
k2 = k2 + 1;
EX(k2) = (en-J+1);
end end end for J = 1 : size(NJ(EX),1)
fprintf(' Столбец %d линейно зависит от столбца %d\n', NJ(J,2), NJ(J,1));
end
X0 = X(:,1:rank(X));
X1 = X(:, rank(X)+1:end);
S = X'*X; %% информационная матрица
% Формирование g-обратной матрицы
A11 = inv(X0'*X0);
A12 = zeros(rank(X), p-rank(X));
A21 = zeros(p-rank(X), rank(X));
A22 = zeros(p-rank(X),p -rank(X) );
Sg = [A11, A12;A21,A22]; % g-обратная матрица
% Определение решения нормального уравнения
H = Sg*S;
z = randn(size(X,2), 1); %%случайный вещественный вектор
P11 = inv(X0’*X0)*X0’*Y;
P21 = zeros(p - rank(X), 1);
% Z11 = zeros(rank(X));
% Z12 = inv(X0’*X0)*X0’*X1;
% Z21 = zeros(p - rank(X), rank(X));
% Z22 = -eye((p - rank(X)));
% A = [P11;P21]+ [Z11,Z12;Z21,Z22]*z
Ep = eye(p);
A = [P11;P21] + [H - Ep]*z;
disp(‘ ‘)
disp(' Оценка вектора параметров A: ')
disp(A)
Результат выполнения программы
Размерность вектора параметров модели: p = 3
Модель наблюдений неполного ранга: rank(X) = 2
Столбец 3 линейно зависит от столбца 1
Оценка вектора параметров A:
54.4373 2.4329
-2.1269
Задание 2 1. Проведите сравнение заданного вектора наблюдений и расчетного вектора наблюдений. Рассчитайте относительную погрешность по каждому значению векторов и определите усредненную
относительную погрешность.
2. Рассчитайте остаточную сумму квадратов RSS.
3. Произведите "зашумление" расчетного вектора наблюдений случайным вектором, распределенным по нормальному закону с параметрами
,
, где — номер компьютера, за которым выполняется лабораторная работа (1, 2, 3, ...). Выполните два предыдущих пункта задания.
4. Выполните решение примера с помощью псевдообратной матрицы
Мура–Пенроуза. Сравните результаты для случая применения - обратной матрицы.
Контрольные вопросы
1. Каковы условия вырожденности информационной матрицы в нормальном уравнении?
2. Что такое остаточная сумма квадратов?
3. Каким свойствам должна удовлетворять псевдообратная матрица
Мура–Пенроуза?
4. Каким основным свойством обладает обобщенная обратная матрица?
5. В каком случае модель будет моделью неполного ранга?
6. Что такое регрессоры?
Регрессионная идентификация линейных непрерывных систем управления
Цель работы: изучить способы оценки (идентификации) параметров
(матриц) непрерывных систем управления на основе регрессионного подхода. Среда программирования — MATLAB.
Теоретическая часть
Описание непрерывных стационарных систем управления в
пространстве состояний имеет вид
(12.1)
где:
— -мерный вектор состояния;
— -мерный вектор управления (входные управляющие воздействия);
—
-мерный вектор выхода системы, – матрица состояния размера
;
— матрица входа размера
;
— матрица выхода размера
;
— матрица обхода размера
Для строго реализуемых систем матрица
В задачу регрессионной оценки (идентификации) линейной стационарной системы (12.1) может входить определение ее структуры и параметров по наблюдаемым данным — входным и выходным сигналам функционирующей системы. Если структура системы задана,
т. е. известны дифференциальные уравнения, описывающие систему, то
в задачу входит определение ее параметров — коэффициентов дифференциальных уравнений — матриц действительных чисел
Для применения регрессионного анализа систему с непрерывным временем (12.1) следует представить в дискретной форме:
(12.2)
где:
— шаг квантования (период дискретизации) по времени;
— целые числа,
,
;
— матрицы дискретной системы тех же размеров, что и для исходной непрерывной системы
[7]
Матрицы
, имеют следующий вид:
(12.3)
(12.4)
Если матрица непрерывной системы не вырожденная, то матрицу можно представить в виде
(12.5)
где — единичная матрица -го порядка
[3]
Поскольку в уравнениях дискретной системы шаг квантования по времени входит слева и справа, то его часто опускают и систему записывают в виде
(12.6)
По известной матрице дискретной системы можно определить матрицу непрерывной системы, логарифмируя обе части уравнения
(12.3). При этом следует иметь в виду, что в (12.3) следует применять матричный экспоненциал, а для обратного преобразования —
матричный логарифм. В системе MATLAB имеются функции матричного экспоненциала и матричного логарифма (см. и
).
В случае неособенной матрицы непрерывной системы уравнение
(12.5) можно разрешить относительно матрицы в виде
(12.7)
Для достаточно малых значений шага квантования можно воспользоваться следующей приближенной формулой:
Таким образом, оценка матриц и непрерывной системы (12.1)
осуществляется через оценку соответствующих матриц дискретной системы (12.2). При этом существуют более универсальные методы
преобразования дискретной системы к непрерывной [
11
,
14
].
Систему разностных уравнений запишем в виде скалярных уравнений:
В задачу входит оценка (идентификация) параметров
,
,
,
Выпишем из (12.8) промежуточное уравнение:
В ходе эксперимента нужно запомнить решений
,
, для идентификации параметров дискретной системы, причем
[1]
. Решения системы (12.8) определяются в результате подачи на нее некоторых входных сигналов (управляющие воздействия). Чтобы проверить, насколько построенная модель точно имитирует или предсказывает данные наблюдений, необходимо сравнить их при одинаковых воздействиях. Эта процедура называется
верификацией модели
[7]
Составим матрицу из элементов размерностью
. Эта матрица будет представлять собой совокупность входных воздействий и выхода системы на момент дискретного времени :
(12.10)
Составим матрицу из искомых коэффициентов уравнения (12.9):
(12.11)
С учетом (12.10), (12.11) запишем уравнение (12.9) в матричном виде:
(12.12)
Сформируем вектор из элементов правой части уравнения (12.12)
после испытаний (значений, решений, наблюдений):
(12.13)
Выражая через
,
, получим
, где
Размерность матрицы равна
Считая, что вектор и матрица известны (входные и выходные сигналы), можно применить метод наименьших квадратов, в соответствии с которым получим следующее нормальное уравнение относительно искомых параметров уравнения (12.9):
(12.15)
Если матрица невырожденная, то оптимальная оценка параметров дискретной системы определяется в виде
(12.16)
Зная
, для всех индексов можно найти коэффициенты
системы уравнений (12.8), т. е. матрицы
,
, а затем, например, по формулам (12.3), (12.7) найти матрицы непрерывной системы
(12.1).
В случае, когда матрица вырожденная или плохо обусловленная,
при решении нормального уравнения (12.6) прибегают к псевдообращению, например, используют псевдообратную матрицу
Мура–Пенроуза.
Расчет вектора состояния дискретной системы можно произвести по соотношению
(12.17)
где
— начальный вектор состояния системы в момент времени,
равный нулю
[3]
Для регрессионной оценки матриц уравнения выхода и используются те же способы и приемы, которые были описаны для оценки матриц и . В случае невырожденной матрицы одновременная оценка матриц и может быть получена с помощью матричного уравнения следующего вида:
(12.18)
где
Примечание. При верификации модели на вход модели и на реальную систему подают сигналы, которые не использовались при идентификации параметров и сравнении выходных сигналов модели и реальной системы.
Практическая часть
В практической части лабораторной работы будут использовать модели стационарных систем управления. Сигналы входа и выхода будут определяться для известной системы, при этом будем считать, что матрицы , уравнения (12.1) и матрицы
, уравнения (12.2)
известны. Но для оценки этих матриц будут использоваться только входные воздействия и значения векторов состояния.
Пример 1. Произведите регрессионную оценку матриц
, системы управления, состоящую из последовательного соединения трех инерционных звеньев с параметрами
. Входное воздействие на систему примите в виде
. Начальное состояние системы примите нулевым по всем переменным состояния.
Изобразим схему объекта управления, представленного через
передаточные функции. Схема показана на рис. 12.1
Рис. 12.1. Схема объекта управления
Связь между входом и выходом, например, первого звена определяется соотношением, выраженным через изображение по Лапласу:
Передаточные функции определены для нулевых начальных условий.
Поэтому комплексную переменную s формально можно заменить на производную, т. е.
Аналогичные преобразования можно провести для других звеньев. В
итоге получим следующую систему дифференциальных уравнений:
(12.21)
Систему (12.21) можно представить в матричном виде. Тогда матрицы системы управления будут иметь следующий вид:
(12.22)
С учетом числовых значений параметров инерционных звеньев будем иметь
(12.23)
Программный код решения примера:
function LAB12;
clc,close all
%% Параметры инерционных звеньев k1 = 1; T1 = 1;
k2 = 2; T2 = 2;
k3 = 3; T3 = 0.5;
%%%------------------------------
%% Матрицы непрерывной системы global A B
A = [-1/T1, 0, 0;
k2/T2, -1/T2, 0;
0, k3/T3, -1/T3]
B = [k1/T1; 0; 0]
C = [0, 0, 1]
D = 0
%%%------------------------------
%% Размерность системы управления n = length(A);
%% Размерность управления r = size(B, 2);
%% Преобразование к дискретной системе
Ts = 0.01; % шаг квантования
%% Расчет матрицы Ad
Ad = expm(A*Ts)
%% Расчет матрицы Bd if abs(det(A)) < 1e-10
Bd = inv(A)*(Ad - eye(size(A)))*B
else syms tau
Bd2 = int(expm(-A*tau)*B, tau, 0, Ts);
Bd = double(Bd2)
end
%%%------------------------------
%% Время наблюдений дискретной системы
Kend = 100; %% число шагов квантования
%% Вид входного воздействия
% % U(t) = exp(-0.5*t).*cos(2*t);
%% Решение разностного уравнения
X0 = zeros(length(Ad),1);
Xk = zeros(length(Ad),Kend);
for k = 1 : Kend sm = zeros(length(Ad), 1);
for J = 0 : k-1
sm = sm + Ad^(k-1-J)*Bd*exp(-0.5*J*Ts)*cos(2*J*Ts);
end
Xk(:,k) = (Ad^k)*X0 + sm;
Uk(1:r,k) = exp(-0.5*k*Ts)*cos(2*k*Ts);
end
%% Хранение массива KSI
KSI = [Xk(:, 2 : end)]';
%% Хранение массива Wk
Wk = [(Xk(:, 1 : end-1))',(Uk(1 : end-1))'];
%% Регрессионная оценка параметров дискретной системы if abs(det(Wk'*Wk)) < 1e-10
F = (pinv(Wk'*Wk)*Wk'*KSI)'; else
F = (inv(Wk'*Wk)*Wk'*KSI)' ; end
%% Идентифицированные матрицы
Adr = F(:, 1 : n)
Bdr = F(:, n+1 : end)
%% Обратное преобразование - оценка матриц А, В
global Areg Breg
Areg = logm(Adr)/Ts if Ts <= 1e-1
Breg = Bdr/Ts elseif Ts > 1e-1 & abs(det(Areg)) <= 1e-10
Breg = inv(inv(Areg)*(Adr - eye(size(Adr))))*Bdr else sd2 = ss(Adr, Bdr, C, D, Ts);
if abs(min(eig(Adr)) ) < 1e-3
sreg = d2c(sd2,'tustin');
else sreg = d2c(sd2,'zox');
end
[Areg, Breg, C, D] = ssdata(sreg) end
%%% Верификация модели
%%% Анализ систем при ступенчатом воздействии global Um
Um = 12;
Xreg0 = zeros(length(Areg), 1);
X0 = zeros(length(A), 1);
T = [0, 12];
[treg, Xreg] = ode23(@fun, T, Xreg0);
[t, X] = ode23(@fun0, T, X0);
h12 = figure(1);
set(h12, 'name','Верификация при ступенчатом воздействии');
line(treg, Xreg, 'linew', 2);
line(t, X,'marker', 'o');
grid on
N = 2*length(Areg);
MN = [ [1:N/2],[1:N/2]]';
str1 = '\bf\it\fontname{times} x\rm\bf_';
str2 = num2str(MN);
str3 = '(\itt\rm\bf)';
legend(strcat(str1, str2, str3), 'location','best');
title('\bf\fontsize{10}\fontname{times}Результат верификации двух моделей систем управления');
xlabel('\bf\fontsize{12}\fontname{times}\it - - - - - - - t - - - - - - -');
ylabel('\bf\fontsize{12}\fontname{times}\it X\rm\bf(\itt\rm\bf) ');
%%%------------------------------ function f = fun0(t,X)
global A B Um f = A*X + B*Um;
function f = fun(tr, Xr)
global Areg Breg Um
f = Areg*Xr + Breg*Um;
В программе преобразование матриц дискретной системы к непрерывной системе может выполняться при ряде условий, в частности, с помощью специализированных функций
(создание класса объекта – дискретной модели),
(преобразование дискретной
модели к непрерывной) и
(извлечение матриц системы). В
функции используются методы преобразования
– экстаполятор нулевого порядка,
– экстраполятор на основе билинейной аппроксимации Тастина [
11
,
14
]. Метод применяется, когда
собственные числа матрицы состояния дискретной системы не лежат вблизи нуля комплексной плоскости.
Результат выполнения программы
A =
-1.0000 0 0 1.0000 -0.5000 0 0 6.0000 -2.0000
B =
1 0
0
C =
0 0 1
D =
0
Ad =
0.9900 0 0 0.0099 0.9950 0 0.0003 0.0593 0.9802
Bd =
0.0101
-0.0001 0.0000
Adr =
0.9900 0.0000 -0.0000 0.0099 0.9950 0.0000 0.0003 0.0593 0.9802
Bdr =
0.0101
-0.0001 0.0000
Areg =
-1.0000 0.0000 -0.0000 1.0000 -0.5000 0.0000 0.0000 6.0000 -2.0000
Breg =
1.0050
-0.0050 0.0001
Как видно из полученных результатов, оценка матриц непрерывной системы является достаточно точной. Верификация модели
производилась при ступенчатом входном воздействии.
Диаграмма переходных процессов показана на рис. 12.2
Рис. 12.2. Переходные процессы по переменным состояния Задание 1 1. Оценку матриц состояния выполните при начальных условиях, равных Х по всем переменным состояния, где — номер компьютера, за которым выполняется лабораторная работа (1, 2, 3, ...). 2. Вычислите регрессионной оценки матриц дискретной системы. 3. Определите приведенную погрешность по компонентам матриц состояния и входа. 4. Постройте усредненную приведенную погрешность по компонентам матриц состояния и входа в зависимости от шага
квантования. Интервал изменения шага квантования примите от до 1.
5. Для оценки матриц непрерывной системы рассмотрите объект управления в виде последовательного соединения двух колебательных звеньев с параметрами
,
,
,
,
,
. Входное воздействие,
используемое для оценки параметров системы, примите в виде
. Верификацию модели произведите при ступенчатом воздействии и нулевых начальных условиях, где —
номер компьютера, за которым выполнятся лабораторная работа (1,
2, 3, ...).
Пример 2. Произведите регрессионную оценку матрицы выхода системы управления, состоящую из последовательного соединения трех инерционных звеньев с параметрами
. Входное воздействие на систему примите в виде
. Начальное состояние системы примите нулевым по всем переменным состояния.
Матрицы системы управления определяются выражениями (12.22).
Проведем регрессионную идентификацию матрицы выхода на основе уравнения
Программный код решения примера:
function LAB122;
clc,close all
%% Параметры инерционных звеньев k1 = 1; T1 = 1;
k2 = 2; T2 = 2;
k3 = 3; T3 = 0.5;
%%%------------------------------
%% Матрицы непрерывной системы global A B
A = [-1/T1, 0, 0;
k2/T2, -1/T2, 0;
0, k3/T3, -1/T3];
B = [k1/T1; 0; 0];
C = [0, 0, 1]
D = 0;
%%%------------------------------
%% Размерность системы управления n = length(A);
%% Размерность управления r = size(B, 2);
%% Размерность выхода m = size(C, 1);
%% Преобразование к дискретной системе
Ts = 0.01; % шаг квантования
%% Расчет матрицы Ad
Ad = expm(A*Ts);
%% Расчет матрицы Bd if abs(det(A)) < 1e-10
Bd = inv(A)*(Ad - eye(size(A)))*B;
else syms tau
Bd2 = int(expm(-A*tau)*B, tau, 0, Ts);
Bd = double(Bd2);
end
%%%------------------------------
%% Время наблюдений дискретной системы
Kend = 100; %% число шагов квантования
%% Вид входного воздействия
% % U(t) = exp(-0.5*t).*cos(2*t);
%% Решение разностного уравнения
X0 = zeros(length(Ad),1);
Xk = zeros(length(Ad),Kend);
for k = 1 : Kend sm = zeros(length(Ad), 1);
for J = 0 : k-1
sm = sm + Ad^(k-1-J)*Bd*exp(-0.5*J*Ts)*cos(2*J*Ts);
end
Xk(:,k) = (Ad^k)*X0 + sm;
Uk(1:r,k) = exp(-0.5*k*Ts)*cos(2*k*Ts);
end
%% Массив выходных значений системы
Ykd = C*Xk; %% как результат измерения
%% Хранение массива KSI
KSI = [Xk(:, 2 : end)]';
%% Хранение массива Wk
Wk = [(Xk(:, 1 : end-1))',(Uk(1 : end-1))'];
%% Регрессионная оценка параметров дискретной системы if abs(det(Wk'*Wk)) < 1e-10
F = (pinv(Wk'*Wk)*Wk'*KSI)'; else
F = (inv(Wk'*Wk)*Wk'*KSI)' ; end
%% Идентифицированные матрицы
Adr = F(:, 1 : n);
Bdr = F(:, n+1 : end);
%% Хранение массива Yk
Yk = Ykd(:, 2:end);
%% Хранение массива YWk
YWk = Xk(:, 1:end-1);
%% Регрессионная оценка матриц дискретной системы if abs(det(YWk*YWk')) < 1e-10
FY = (pinv(YWk*YWk')*YWk*Yk'); else
FY = (inv(YWk*YWk')*YWk*Yk'); end
%% Идентифицированная матрица выхода
Cdr = FY' %% для дискретной системы
Creg = Cdr %% для непрерывной системы
%% Обратное преобразование - оценка матриц А, В
global Areg Breg
Areg = logm(Adr)/Ts;
if Ts <= 1e-1
Breg = Bdr/Ts;
elseif Ts > 1e-1 & abs(det(Areg)) <= 1e-10
Breg = inv(inv(Areg)*(Adr - eye(size(Adr))))*Bdr; else sd2 = ss(Adr, Bdr, C, D, Ts);
if abs(min(eig(Adr)) ) < 1e-3
sreg = d2c(sd2,'tustin');
else sreg = d2c(sd2,'zox'); end
[Areg, Breg, C, D] = ssdata(sreg); end
%%% Верификация модели
%%% Анализ систем при ступенчатом воздействии
global Um
Um = 12;
Xreg0 = zeros(length(Areg), 1);
X0 = zeros(length(A), 1);
T = [0, 12];
[treg, Xreg] = ode23(@fun, T, Xreg0);
Yreg = Xreg*Creg'; %% выход системы
[t, X] = ode23(@fun0, T, X0);
Y = X*C'; %% выход системы h12 = figure(1);
set(h12, 'name','Верификация при ступенчатом воздействии');
line(treg, Yreg, 'linew', 2);
line(t, Y, 'marker', 'o','color','r');
str1 = ...
'\bf\fontsize{11}\fontname{times}\itY_r_e_g\rm\bf(\itt\rm\bf)';
str2 ='\bf\fontsize{11}\fontname{times}\itY\rm\bf(\itt\rm\bf)';
legend(str1, str2, 'location','best');
title('\bf\fontsize{10}\fontname{times}Результат верификации двух моделей по выходу');
xlabel('\bf\fontsize{12}\fontname{times}\it - - - - - - - t - - - - - - -');
ylabel('\bf\fontsize{12}\fontname{times}\it Y\rm\bf(\itt\rm\bf) '); grid on
%%%------------------------------ function f = fun0(t,X)
%% Для заданной системы global A B Um f = A*X + B*Um;
function f = fun(tr, Xr)
%% Для идентифицированной системы global Areg Breg Um f = Areg*Xr + Breg*Um;
Результат выполнения программы C = 0 0 1 Cdr = 0.0003 0.0592 0.9802 Creg = 0.0003 0.0592 0.9802 Диаграмма выходного процесса системы при верификации показана на рис. 12.3 Рис. 12.3. Переходный процесс по выходу систем
Задание 2 1. Постройте усредненную погрешность оценки выходной матрицы при изменении шага квантования от до 1.
2. Произведите регрессионную оценку выходной матрицы при размерности выхода, равной двум (
). Постройте также переходные процессы по выходным переменным.
Контрольные вопросы
1. Что такое матричный экспоненциал?
2. Какие виды погрешностей обусловливают оценку матриц линейной модели системы управления?
3. Что такое верификация?
4. Как называется уравнение, которое используется для регрессионной оценки матриц линейной системы управления?
5. В каком случае и на каком этапе используется псевдообратная матрица Мура–Пенроуза при регрессионной идентификации матриц линейной системы управления?
Список литературы
1. Асатурян В.И. , Теория планирования эксперимента, М.: Радио и связь, 1983. 248 с
2. Афонин В.В., Мурюмин С.М., Федосин С.А, Основы анализа систем массового
обслуживания, Саранск: Изд-во Мордов. ун-та, 2003. 236 с
3. Афонин В.В., Федосин С.А., Иконников С.Е, Основы теории управления. Лабораторный
практикум, Саранск: Изд-во Мордов. ун-та, 2008. 244 с
4. Вентцель Е.С., Овчаров Л.А, Теория случайных процессов и ее инженерные приложения, М.:
Высшая школа, 2000. 383 с. 2-е изд
5. Гаврилов В.Р., Иванова Е.Е., Морозова В.Д, Кратные и криволинейные интегралы.
Элементы теории поля, М.: Изд-во МГТУ им. Н.Э. Баумана, 2003. 496 с
6. Гмурман В.Е, Теория вероятностей и математическая статистика, М.: Высшая школа,
2003. 479 с. 9-е изд
7. Гроп Д, Методы идентификации систем, М.: МИР, 1979. 302 с
8. Иглин С.П, Математические расчеты на базе MATLAB, СПб.: БХВ–Петербург, 2005. 640 с
9. Кнут Д, Искусство программирования для ЭВМ. Т. 2: Получисленные алгоритмы, М.:
Издательский дом "Вильямс", 2003. 832 с. Изд. 3 10. Кудрявцев Л.Д, Курс математического анализа (в 2-х томах), М.: Высшая школа, 1981, Т.
I. 687 с
11. Лазарев Ю, Моделирование процессов и систем в MATLAB, СПб.: Питер; Киев:
Издательская группа BXV, 2005. 512 с
12. Львович Я.Е., Фролов В.Н, Теоретические основы конструирования, технологии и
надежности РЭА, М.: Радио и связь, 1986. 192 с
13. Горяинов В.Б., Павлов И.В., Цветкова Г.М. и др, Математическая статистика, М.: Изд- во МГТУ им. Н.Э. Баумана, 2001. 424 с
14. Медведев В.С., Потемкин В.Г, Control System Toolbox. MATLAB 5 для студентов, М.:
ДИАЛОГ – МИФИ, 1999. 287 с
15. Мэтьюз Д.Г, Численные методы. Использование MATLAB, Издательский дом "Вильямс",
2001. 720 с. 3-е издание
16. Рыжиков Ю.И, Имитационное моделирование. Теория и технологии, СПб.: КОРОНА
принт; М.: Альтекс-А, 2004. 384 с
17. Себер Дж, Линейный регрессионный анализ, М.: Мир, 1980. 456 с
18. Сигорский В.П, Математический аппарат инженера, Киев: Технiка, 1975. 768 с
19. Советов Б.Я., Яковлев С.А, Моделирование систем: Учеб. для вузов, М.: Высшая школа,
2005. 343 с. 4-е изд., стереотип
20. Таха Х, Введение в исследование операций, М.: Мир, 1985. 496 с
21. Томашевский В., Жданова Е, Имитационное моделирование в среде GPSS, М.: Бестселлер,
2003. 416 с
|
Скачать 5.72 Mb.