Рис. 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;
