(2.9)

где

Таким образом, получены две эквивалентные формы описания системы: дифференциальное уравнение (2.2) с начальными условиями (2.3) и интегральное уравнение (2.9). Функция  в выражении (2.9) представляет собой полином, коэффициенты которого зависят от начальных условий (2.3) и от множества  искомых параметров настройки системы автоматического управления (регулирования). Перепишем , изменив порядок суммирования

Введем следующие обозначения:

Тогда полином  можно записать следующим образом

где - вектор-столбец начальных условий; - вектор-столбец полиномов .

Рассмотрим левую часть уравнения (2.9). Представим функции, входящие в нее, в виде разложений в ряд по ортонормированному базису .

Имеем

, (2.10)

где  - спектральная характеристика выходного сигнала , элементы которой определяются из соотношения

 (2.11)

где  - квадратная матрица размерностью , элементы которой определяются из выражения

Подставив полученные разложения (2.10) и (2.11) в левую часть уравнения (2.9) и учитывая, что , где - единичная, в силу ортонормированности базисных функций, получим

 (2.12)

где  - матрица спектральной характеристики инерционной части системы размерностью .

Сделаем аналогичные преобразования для правой части уравнения (2.9).

, (2.13)

где  - спектральная характеристика сигнала на входе системы, элементы которой определяются из соотношения

 (2.14)

где  - квадратная матрица размерностью  спектральной характеристики форсирующей части системы, элементы которой определяются из выражения

 (2.15)

где  - матрица размерностью  элементы которой определяются из соотношения

Подставляя разложения (2.13), (2.14) и (2.15) в (2.9) и делая соответствующие преобразования, получим

 (2.16)

Таким образом, уравнение (2.9) с учетом (2.12) и (2.16) можно переписать в следующем виде

 (2.17)

Рассмотрим теперь функционал (2.4). Имеем

Так как , то последние выражение можно записать в следующем виде

 (2.18)

или

где

. (2.19)

Здесь спектральная характеристика эталонного сигнала  или задана или, в случае задании эталонного сигнала , определяется из выражения

, .

Таким образом, задача определения входного сигнала  (точнее множества ) и множества  неизвестных параметров настройки системы управления (2.2), (2.3) сводиться к задаче безусловной минимизации функционала (2.18) по элементам множеств  и , т.е.

.

На рисунке 2.1 представлена структурная схема алгоритма решения поставленной задачи.

Рис 2.1 Структурная схема алгоритма решения обратной задачи динамики спектральным методом

4. Практическая часть

Рассмотрим отдельный блок системы самонаведения, структурная схема которого представлена на рисунке 1.

Рис. 1. Структурная схема системы

Задан эталонный закон изменения угла , график которого представлен на рисунке 2.

Рис. 2. График эталонного закона изменения угла

Задача формулируется следующим образом. Необходимо найти управление  такое, которое обеспечит на выходе сигнал , максимально близкий к заданному эталонному закону.

5. Практическая часть

Данная задача относится к разряду неккоректных и мы будем решать её с применением оптимизационных методов.

Для решения данной задачи воспользуемся методом матричных операторов. В этом случае структурную схему можно представить в следующем виде (рис. 3).

Рис. 3. Структурная схема системы в операторной форме

В качестве ортонормированной системы использовалась система функций Уолша с удержанием  элементов. В этом случае матричные операторы основных элементов системы будут следующими (представлены подматрицы размерностью ):

;

;

;

.

Спектральная характеристика сигнала  следующая (представлены первые пять элементов):

.

Решение поставленной задачи будем выполнять в следующие два этапа.

1. Поскольку известен эталонный выходной сигнал, то из уравнения

можно найти спектральную характеристику эталонного сигнала на выходе нелинейного элемента. Решая уравнение относительно коэффициентов  с использованием метода Гаусса-Ньютона получены следующие числовые значения коэффициентов:

.   

График соответствующего сигнала представлен на рисунке 4.

Рис. 4. График сигнала, который необходимо получить на выходе нелинейного элемента

Однако на выходе нелинейного элемента можно получить сигнал, представленный на рисунке 5 (ниже показаны первые пять элементов спектральной характеристики).

Рис. 5. Реальный сигнал на выходе нелинейного элемента

.

Тогда из находим эталонный сигнал на выходе, который может обеспечить данная система (рис. 6). Его спектральная характеристика:

.      

Рис. 6. Графики требуемого эталонного сигнала и эталонного сигнала, который можно получить

2. В результате решения предыдущего этапа найдены спектральные характеристики эталонного выходного сигнала, который может обеспечить данная система, и эталонного сигнала, которой необходимо получить на входе нелинейного элемента.

Далее искомый сигнал  представим в виде

,  

где  некоторая система линейно независимых функций.

В результате можно для спектральной характеристики сигнала на входе нелинейного элемента записать следующую зависимость.

,

где  – спектральная характеристика -го элемента системы . Поскольку известны спектральные характеристики эталонных сигналов  и , то между левой и правой частями выражения будет иметь место невязка

,       

зависящая от неизвестных коэффициентов , . Сформировав функционал

,    

исходную задачу синтеза входного сигнала можно свести к задаче поиска минимума функционала на множестве допустимых значений коэффициентов , , т.е.

.

При решении задачи в качестве системы функций  использовались экспоненциальные функции: . Минимум функционала искался с использование алгоритма Нелдера-Мида (алгоритма безусловной минимизации). В качестве начальных значений искомых коэффициентов были приняты нулевые. При этом значение функционала :

.

Были получены следующие оптимальные значения искомых коэффициентов:

;

;

;

;

;

;

;

;

;

.

Значение функционала в оптимальной точке:

.

Следовательно, входной сигнал имеет следующий вид:

.

На рисунке 7 представлен график сигнала .

Рис. 7. График синтезируемого входного сигнала

На рисунке 8 представлены результаты анализа системы с использованием метода Рунге-Кутта для найденного входного сигнала и для сравнения приведены графики требуемого эталонного выходного сигнала и эталонного сигнала, который может обеспечить данная система.

Рис. 8. Графики выходных сигналов системы

Таким образом, можно построить следующий алгоритм решения задачи синтеза входного сигнала нелинейной системы:

1) задается эталонный выходной сигнал;

2) из находится сигнал на выходе нелинейного элемента, который на выходе системы обеспечивает требуемый эталонный процесс;

3) найденный в предыдущем пункте сигнал представляется как сигнал на входе нелинейного элемента и находится реальный сигнал на выходе нелинейного элемента и уточняется эталонный сигнал на выходе системы;

4) поскольку известны сигналы на входе нелинейного элемента и на выходе системы, то, представив искомый входной сигнал в виде , строится невязка и функционал ;

5) минимизируя полученный функционал, находятся числовые значения искомых коэффициентов , ;

6) проводится анализ полученных результатов.

5. Результаты расчёта

1. Эталонный закон изменения угла teta(t)

Число точек квантования по времени: Nt = 499;

Шаг квантования: h_t = 0.020000 c;

Время поражения цели: T = 9.960000 c;

2. Числовые значения параметров системы самонаведения

Krp = 1.000000;

Trp = 0.330000, с;

Xmax = 0.418879, рад;

Ksn = 0.283000, рад/с;

Tsn = 0.155000, с;

DZsn = 0.052000;

V = 686.700000, м/с;

G = 9.810000, м/с^2;

Kdy = 0.140000;

Kv = 1.200000, c;

mu = 0.115000, с;

Tc = 3.050000, с;

3. Базис - функции Уолша

Число элементов: Nl = 64;

Оператор интегрирования Ai размерностью 5x5

+4.980000e+000 +2.490000e+000 +0.000000e+000 +1.245000e+000 +0.000000e+000

-2.490000e+000 +0.000000e+000 +1.245000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 -1.245000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

-1.245000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +6.225000e-001

+0.000000e+000 +0.000000e+000 +0.000000e+000 -6.225000e-001 +0.000000e+000

Оператор дифференцирования Ad размерностью 5x5

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

4. Матричные операторы системы

Arp размерностью 5x5

+9.668675e-001 +3.313252e-002 -3.310223e-002 +3.310224e-002 -3.171612e-002

-3.313252e-002 +9.006024e-001 +9.930669e-002 +3.310223e-002 -3.171610e-002

-3.310223e-002 -9.930669e-002 +8.343980e-001 +3.307196e-002 -3.168711e-002

-3.310224e-002 +3.310223e-002 -3.307196e-002 +7.682541e-001 +2.220418e-001

-3.171612e-002 +3.171610e-002 -3.168711e-002 -2.220418e-001 +7.048218e-001

Asn размерностью 5x5

+2.824568e-001 -1.904545e-003 -3.561384e-003 +6.907620e-003 -4.945520e-003

+1.904545e-003 +2.862659e-001 +1.403039e-002 +3.561384e-003 -6.087807e-003

-3.561384e-003 -1.403039e-002 +2.791431e-001 +1.571978e-002 -7.488732e-003

-6.907620e-003 +3.561384e-003 -1.571978e-002 +2.477036e-001 +3.501836e-002

-4.945520e-003 +6.087807e-003 -7.488732e-003 -3.501836e-002 +2.378125e-001

Aos1 размерностью 5x5

-6.831527e+001 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 -6.831527e+001 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 -6.831527e+001 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 -6.831527e+001 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 -6.831527e+001

Aos2 размерностью 5x5

+9.800000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +9.800000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +9.800000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +9.800000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 +9.800000e+000

Apr размерностью 5x5

+2.730338e-001 +9.500096e-003 -1.194782e-002 +1.128111e-002 -7.250387e-003

-9.500096e-003 +2.540337e-001 +3.517674e-002 +1.194782e-002 -8.567413e-003

-1.194782e-002 -3.517674e-002 +2.301380e-001 +1.306212e-002 -5.530241e-003

-1.128111e-002 +1.194782e-002 -1.306212e-002 +2.040138e-001 +5.461586e-002

-7.250387e-003 +8.567413e-003 -5.530241e-003 -5.461586e-002 +1.895130e-001

Aos размерностью 5x5

-5.851527e+001 +0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 -5.851527e+001 +0.000000e+000 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 -5.851527e+001 +0.000000e+000 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 -5.851527e+001 +0.000000e+000

+0.000000e+000 +0.000000e+000 +0.000000e+000 +0.000000e+000 -5.851527e+001

As размерностью 5x5

+3.194591e-001 +1.707523e-001 +5.751752e-004 +8.508857e-002 +5.722004e-004

-1.707523e-001 -2.204553e-002 +8.393822e-002 -5.751752e-004 +5.722004e-004

+5.751752e-004 -8.393822e-002 -2.089518e-002 -5.751662e-004 +5.721915e-004

-8.508857e-002 -5.751752e-004 +5.751662e-004 -1.974485e-002 +3.882646e-002

+5.722004e-004 -5.722004e-004 +5.721915e-004 -3.882646e-002 -1.860045e-002

5. СХ эталонного выхода

Ctheta размерностью 5x1

+2.948462e-001

-7.002572e-002

-4.945100e-002

-5.104576e-002

-1.450117e-002

Начальные значения искомых коэффициентов

Cu_0 размерностью 5x1

+0.000000e+000

+0.000000e+000

+0.000000e+000

+0.000000e+000

+0.000000e+000

Oshibka_0 = 3.145671e-001

Conditioning of Gradient Poor - Switching To LM method

Optimization terminated: directional derivative along

search direction less than TolFun and infinity-norm of

gradient less than 10*(TolFun+TolX).

Оптимальные значения искомых коэффициентов

Cu_opt размерностью 5x1

+5.349004e-001

+5.156158e-001

+3.167675e-001

+3.345843e-001

+3.459092e-002

Oshibka_0 = 1.444098e-004

-------------------------------------------------------------

Время расчета:

0 часов, 0 минут, 34.703 секунд.

Приложение

1) % Программа синтеза управления системы самонаведения (рассматривается часть % системы) методами обратных задач динамики с использованием метода % матричных операторов (линейная модель)

close all;

clear all;

clc;

my_tic;

global Nl;

global U tgl;

global Krp Trp Ksn Tsn DZsn V G Kdy Kv mu Tc Xmax;

%% 1. Эталонный закон изменения угла teta(t)

% Время наведения

fId = fopen('t_navedenija.dat','r');

t_f = fread(fId,inf,'real*8')';

fclose(fId);

Nt_f = length(t_f);

h_t_f = t_f(2)-t_f(1);

T = t_f(Nt_f);

% угол theta(t)

fId = fopen('theta_navedenija.dat','r');

theta_f = fread(fId,[1 Nt_f],'real*8');

fclose(fId);

% расстояние до цели

fId = fopen('r_navedenija.dat','r');

r_f = fread(fId,[1 Nt_f],'real*8');

fclose(fId);

fprintf('1. Эталонный закон изменения угла teta(t)\n');

fprintf('Число точек квантования по времени: Nt = %i;\n',Nt_f);

fprintf('Шаг квантования: h_t = %f c;\n',h_t_f);

fprintf('Время поражения цели: T = %f c;\n',T);

fprintf('\n');

my_plot2(t_f,theta_f,'t, c','theta(t), рад');

my_plot2(t_f,r_f,'t, c','r(t), м');

% пересчет на больший шаг квантования

Nt = 64;

h_t = T/(Nt-1);

t = 0: h_t: T;

theta = spline(t_f,theta_f,t);

r = spline(t_f,r_f,t);

my_plot2(t,theta,'t, c','theta(t), рад');

my_plot2(t,r,'t, c','r(t), м');

%% 2. Параметры системы

% Числовые значения параметров системы самонаведения

Krp = 1; %

Trp = 0.33; % с

Xmax = 24*pi/180; % рад

Ksn = 0.283; % рад/с

Tsn = 0.155; % с

DZsn = 0.052; %

V = 70*9.81; % м/с

G = 9.81; % м/с^2

Kdy = 0.14; %

Kv = 1.2; % c

mu = 0.115; % с

Tc = 3.05; % с

fprintf('2. Числовые значения параметров системы самонаведения\n');

fprintf('Krp = %f;\n',Krp);

fprintf('Trp = %f, с;\n',Trp);

fprintf('Xmax = %f, рад;\n',Xmax);

fprintf('Ksn = %f, рад/с;\n',Ksn);

fprintf('Tsn = %f, с;\n',Tsn);

fprintf('DZsn = %f;\n',DZsn);

fprintf('V = %f, м/с;\n',V);

fprintf('G = %f, м/с^2;\n',G);

fprintf('Kdy = %f;\n',Kdy);

fprintf('Kv = %f, c;\n',Kv);

fprintf('mu = %f, с;\n',mu);

fprintf('Tc = %f, с;\n',Tc);

fprintf('\n');

%% 3. Формирование ортонормированного базиса

Nl = Nt;

setsize(Nl);

settime(T);

Ai = mkint; % оператор интегрирования

Ad = inv(Ai); % оператор дифференцирования

Ae = eye(Nl); % единичная матрица

fprintf('3. Базис - функции Уолша\n');

fprintf('Число элементов: Nl = %i;\n',Nl);

pr_matrix(Ai,'Оператор интегрирования Ai')

pr_matrix(Ad,'Оператор дифференцирования Ad')

%% 4. Расчет операторов системы

Arp = inv(Trp*Ae+Ai)*(Krp*Ai);

Asn = inv(Tsn^2*Ae+2*DZsn*Tsn*Ai+Ai*Ai)*(Ksn*Ai*Ai);

Aos1 = Kv*mu*Tc*Ad*Ad+Kv*(mu+Tc)*Ad+Kv*Ae;

Aos2 = (Kdy*V/G)*Ae;

Apr = Asn*Arp;

Aos = Aos1+Aos2;

As = inv(Ae+Aos*Apr)*Apr;

As = Ai*As;

fprintf('4. Матричные операторы системы\n');

pr_matrix(Arp,'Arp');

pr_matrix(Asn,'Asn');

pr_matrix(Aos1,'Aos1');

pr_matrix(Aos2,'Aos2');

pr_matrix(Apr,'Apr');

pr_matrix(Aos,'Aos');

pr_matrix(As,'As');

%% 5. Расчет спектральной характеристики эталонного выхода

Ctheta = fwht(theta');

fprintf('5. СХ эталонного выхода\n');

pr_matrix(Ctheta,'Ctheta');

%% 6. Синтез входного сигнала

Cu_0 = zeros(Nl,1);

fprintf('Начальные значения искомых коэффициентов\n');

pr_matrix(Cu_0,'Cu_0');

oshibka = sqrt((As*Cu_0-Ctheta)'*(As*Cu_0-Ctheta));

fprintf('Oshibka_0 = %e\n',oshibka);

my_function = @(Cu)sqrt((As*Cu-Ctheta)'*(As*Cu-Ctheta));

% optimset('Display','iter','NonlEqnAlgorithm','gn','TolFun',1e-8,...

Cu = fsolve(my_function,Cu_0,...

 optimset('NonlEqnAlgorithm','gn','TolFun',1e-8,...

 'TolX',1e-8,'MaxFunEvals',50000,'MaxIter',50000));

% Cu = inv(As)*Ctheta;

fprintf('Оптимальные значения искомых коэффициентов\n');

pr_matrix(Cu,'Cu_opt');

oshibka = sqrt((As*Cu-Ctheta)'*(As*Cu-Ctheta));

fprintf('Oshibka_0 = %e\n',oshibka);

U = iwht(Cu)';

tgl = t;

my_plot2(t,U,'t, c','U(t)');

%% 7. Анализ полученных результатов (метод Рунге-Кутта (ode45))

[tt,yy] = ode45(@ode_navedenija1,t,[0 0 0 0]);

theta_rr = yy(:,1)';

my_plot2(t,[theta;theta_rr],'t, c','theta(t), рад','',['эталонный ';'реальный ']);

my_toc;

2) второстепенные программы:

function dy = ode_navedenija1(t,y);

global U tgl;

global Krp Trp Ksn Tsn DZsn V G Kdy Kv mu Tc Xmax;

a32 = -1/(Tsn^2);

a33 = -2*DZsn/Tsn;

a3f = Ksn/(Tsn^2);

a42 = -(Krp/Trp)*(Kv-Kv*mu*Tc/(Tsn^2)+Kdy*V/G);

a43 = -(Krp/Trp)*(Kv*(mu+Tc)-2*Kv*mu*Tc*DZsn/Tsn);

a44 = -1/Trp;

a4f = -(Krp/Trp)*Kv*Ksn*mu*Tc/(Tsn^2);

b4 = Krp/Trp;

u = spline(tgl,U,t);

dy = zeros(4,1);

dy(1) = y(2);

dy(2) = y(3);

y4 = y(4);

dy(3) = a32*y(2)+a33*y(3)+a3f*y4;

dy(4) = b4*u+a42*y(2)+a43*y(3)+a44*y(4)+a4f*y4;

Страницы: 1, 2