t-anv

Обратные задачи динамики
Что-то делаю не так. Система сходиться гдето за 1000с, а должна за 30с.
В Матлабе
Метод:
Концепция обратных задач динамики. Алгоритм назначенных траекторий движения

Пусть движение управляемой системы подчиняется дифференциальному уравнению
x ̇(t)=f(x,t)+g(u,t) (3.1.1)
где x=‖x_1…x_n ‖^T – вектор координат состояния; u=‖u_1…u_m ‖^T - вектор управляющих функций; f(x,t),g(u,t) – заданные функции. В пространстве состояний задана некоторая траектория
x^* (t)=ψ(t), 0<t<T, 0<T<∞. (3.1.2)
Путем выбора функций u_R (t) требуется найти такие управляющие силы g_v (u,t), при которых движение системы (3.1.1) точно или приближенно следует по траектории ψ(t).
При условии, что поставленная задача имеет точное решение, искомый вектор управлений определяется из уравнения
g(u,t)=ψ ̇(t)-f(ψ(t),t). (3.1.3)
Найденные отсюда управляющие функции u_1 (t),…,u_m (t) являются программными, так как они изменяются в соответствии с заданной вектор-функцией (3.1.2). Система управления в таком случае является разомкнутой, ее структурная схема состоит из последовательного соединения регулятора и управляемого объекта (рисунок 3.1.1).

Рис.3.1.1. Разомкнутая система управления
В качестве примера рассмотрим простейший случай, когда управляемый процесс подчиняется дифференциальному уравнению
x ̇(t)=ax(t)+bu(t),x(0)=x_0,
где x – скалярная функция. Пусть требуется осуществить движение по назначенной траектории ψ(t), определяемой равенством
ψ(t)=x_0 e^(-ct),c>0.
Чтобы найти управляющую функцию, реализующую движение по назначенной траектории, подставим выражение для ψ(t) в уравнение движения вместо x(t). Имеем
u(t)=-b^(-1) (c+a)x_0 e^(-ct).
Найденная управляющая функция является программной, так как она определена как функция времени. Схема реализации этого управления показана на рисунке 3.1.2.

Рис.3.1.2.
Можно видеть, что для формирования управляющей функции u(t) необходимо генерировать функцию времени exp⁡(-ct), которая определяет назначенную траекторию. С учетом найденного управления процесс в рассматриваемой системе описывается неоднородным уравнением
x ̇(t)=ax(t)-(a+c) x_0 e^(-ct),x(0)=x_0.
Решение этого уравнения x_0 e^(-ct) в точности совпадает с назначенной траекторией ψ(t).
Уравнение (3.1.3) в интересных для практических приложений случаях оказывается неразрешимым относительно управляющих функций u_R (t). Задача обладает лишь кажущейся простотой. Ниже указаны возможные варианты построения программных управлений применительно к математической модели движения управляемого объекта следующего вида:
x(t)=f(x,t)+∑_(R=1)^m▒c_R u_R (t), (3.1.4)
где c_R n – мерные векторы параметров; u_R – скалярные функции – искомые управления.
Определим вектор-функцию
z(t)=x(t)-ψ(t). (3.1.5)
Тогда вместо (3.1.4) можно записать дифференциальное уравнение для новой переменной
z ̇(t)=f(z,t)+∆(t), (3.1.6)
где приняты следующие обозначения:
f(z,t)=f(z+ψ,t)-f(ψ,t), (3.1.7)
∆(t)=f(ψ,t)-ψ ̇(t)+∑_(R=1)^m▒c_R u_R (t).
Из (3.1.4), (3.1.7) заключаем, что задача осуществления заданной траектории сводится к аппроксимации функции
v(t)=ψ ̇(t)-f(ψ,t)
линейной конструкцией
u(t)=∑_(R=1)^m▒c_R u_R (t). (3.1.8)
Действительно, непосредственная подстановка (3.1.5) в уравнение (3.1.4) дает
ψ ̇(t)=f(ψ,t)+∑_(R=1)^m▒c_R u_R (t).
Отсюда следует, что в случае v(t)≡u(t), т.е. при
ψ ̇(t)-f(ψ,t)≡∑_(R=1)^m▒c_R u_R (t), (3.1.9)
Имеет место точная реализация траектории ψ(t)≡x(t).
В практических случаях условие (3.1.9), как правило, недостижимо. Поэтому назначенную траекторию можно осуществить с определенной степенью приближения. Можно приближенно строить вектор управления (3.1.8) из условия реализации минимума величины
δ=(∫_0^T▒〖∆^T (t)∆(t)dt〗)^(1/2), (3.1.10)
где символ «т» означает операцию транспонирования. В такой постановке возможны три варианта построения управления.
Вариант 1. Заданы функции u_R (t). Требуется найти векторы c_R таким образом, чтобы δ, определяемая равенством (3.1.10), была минимальна.
Вариант 2. Заданы векторы c_R. Требуется найти систему функций u_1 (t),…,u_m (t), при которых реализуется минимум (3.1.10).
Вариант 3. Найти оптимальную систему векторов c_R и функций u_R (t), при которых достигается минимум оценки (3.1.10).
Необходимо отметить, что при обсуждении рассматриваемой задачи предполагается, что начальные точки действительной x(t) и назначенной x*(t) траекторий совпадают. Если начальное состояние управляемой системы не соответствует желаемому, то следует предварительно осуществить переходный процесс.
Синтез системы автоматического регулирования заключается в построении обратных связей, реализующих расчетную или априорную программу движения. Под программами понимаются «как запрограммированные динамические характеристики системы или программа ее функционирования, так и изменение выходных координат системы по заданной функциональной зависимости».
Можно видеть, таким образом, что речь идет здесь о синтезе замкнутых автоматических систем с позиции обратных задач динамики: задана программа движения управляемого объекта, требуется синтезировать такую систему, движение которой осуществляется по предписанной программе.
Рассмотрим процедуру построения вычислительного алгоритма дискретного управления.
Задана модель движения системы
x ̇(t)=f(x,u), (3.1.11)
где x, u – скалярные функции. Относительно функции f(x,u) предполагается, что она дифференцируема по u. Задача ставится следующим образом. Пусть в текущий момент времени t известны требуемые значения x*(t+T/2), x*(t+T) регулируемой переменной x, соответствующие моментам времени t+T/2, t+T. Здесь T – достаточно малая величина, так что с необходимой степенью приближения на основание (3.1.11) можно записать
x(t+T)=x(t)+f[t+T/2]T, (3.1.12)
где для краткости обозначено
f[t+T/2]=f[x(t+T/2),u(t+T/2)].
Требуется найти такое приращение ∆u(t) управляющей функции u(t): при котором в моменты времени t+T/2, t+T переменная x примет заданные значения.
Представление управляемого процесса в виде (3.1.12) означает, что управляющая функция вычисляется в дискретные моменты времени t=kT, k=0, 1, …. В интервалах (t,t+T) ее величина остается постоянной, что отражено на рисунке 3.1.3.

Рис.3.1.3. Представление управляемого процесса
Поскольку соотношение (3.1.12) лишь приближенно описывает течение управляемого процесса, то можно говорить о приближенной реализации заданной траектории движения x*(t+T).
Для определения искомой величины ∆u(t) поступаем следующим образом. Согласно (3.1.12) и рисунку 3.1.3 можно записать
x(t+T)=x(t)+f[t-T/2]T+∆f[t] T/2. (3.1.13)
Если f(x,u) есть медленно меняющаяся функция, то величину ∆f[t] можно приближенно определить следующим образом:
∆f[t]≅df(x,u)/du ∆u≅∆f[t-T]/∆u(t-T) u(t).
Далее, в соответствии с (3.1.12) справедливы равенства
f[t-T/2]=2/T [x(t)-x(t-T/2)];
∆f[t-T]=f[t-T/2]-f[t-3/2 T]= (3.1.14)
=2/T [x(t)-2x(t-T/2)+x(t-T)].
Кроме того, аналогично (3.1.14) можно написать
∆f[t]=2/T [x(t)-2x(t+T/2)+x(t+T)]. (3.1.15)
Следовательно, из формулы
∆u(t)=∆u(t-T)(∆f[t])/(∆f(t-T))
С учетом (3.1.14) и (3.1.15) находим приращение управляющей функции
∆u(t)=∆u(t-T) (x(t)-2x^* (t+T/2)+x^* (t+T))/(x(t)-2x(t-T/2)+x(t-T) ). (3.1.16)
Заметим, что это выражение получено при условии, что в (3.1.15) принято
x(t+T)=x^* (t+T),x(t+T/2)=x^* (t+T/2)
Иначе говоря, приращение ∆u(t) найдено из условия, чтобы в моменты времени t+T/2, t+T регулируемая переменная принимала заданные значения, как это требуется по формулировке задачи.
Таким образом, соотношение (3.1.16), дополненное равенством
u(t)=u(t-T)+∆u(t), (3.1.17)
составляет содержание алгоритма формирования управляющей функции. Как следует (3.1.16) и (3.1.17), для вычисления u(t) необходимо знать состояние управляемого объекта для трех моментов времени:
x(t),x(t-T/2),x(t-T).
Кроме того, в памяти вычислителя должны сохраняться значения
u(t-T),∆u(t-T).
Формулы (3.1.16) и (3.1.17) могут быть записаны в рекуррентной форме.
Для этой цели примем
T^'=T/2,t=kT^',k=0,1,… (3.1.18)
В таком случае алгоритм вычисления управляющей функции принимает вид
∆u(kT^' )=∆u[(k-2) T^' ] (x^* [(k+2) T^' ]-〖2x〗^* [(k+1) T^' ]+x^* (kT^' ))/(x[(k-2) T^' ]-2x[(k-1) T^' ]+x(kT^') );
u(kT^' )=u[(k-2) T^' ]+∆u(kT^' ),k=0,1,… (3.1.19)
Согласно (3.1.18) период дискретности определяется величиной T’.
Итак, для нахождения управляющего воздействия будем использовать рекуррентное соотношение (3.1.19) и условие (3.1.18).


Моя система:
W=(2.7488*10^(-4)s^2+2.268*10^(-5)s+1.18*10^(-7))/(s^3+0.0733s^2+0.001475s+5.25*10^(-6))

Эталонный сигнал, за который осуществляется слежение, имеет вид:
Xэталон.вода(t)=30e^(-0.02t)-9e^(-0.033t)-21e^(-0.04t) . Это x*(kT'). Если вместо t подставить kT'

Листинг программы:

clear all
close all
clc
% Обявляем переменные
global a1 a2 a3 F1 F2 F3 un
% Задаем параметры системы
a1=0.0733;
a2=0.001475;
a3=5.25e-6;
F1=2.7488e-4;
F2=2.5313e-6;
F3=-4.7298e-7;
% T1=0.5;
% T2=0.25;
% T3=0.05;
% T4=1;
% K=1;
% Задаем интервал времени (10с)
Tend=500; % Задаем шаг
h=1;
% Задаем Н.У.
x(:,1) = [2; 0; 0];
x(:,2) = [-0.15; 0; 0] ;
x(:,3) = [0; 0; 0] ;
% Задаем начальное значение приращения управляющей функции
deltaU = [0 -1.8];
% Задаем начальное значение управлякщей функции
u = [2 0 0];
% Задаем начальное значение входного воздействия
for i=1:4
y(i)=30*exp(-0.02*i*h)-9*exp(-0.033*i*h)-21*exp(-0.04*i*h);
% y(i) = 5.6*exp(-0.04*i*h)+2.4*exp(-0.033* i*h)-8*exp(-0.02*i*h));
end
% Реализуем сам алгоритм нахождения приращения управлякщей функции и самой управлякщей функции
for k=3:Tend/h
% y(k+2)= 5.6*exp(-0.04*(k+2)*h)+2.4*exp(-0.033*(k+2)*h)-8*exp(-0.02*(k+2)*h);
y(k+2)=30*exp(-0.02*(k+2)*h)-9*exp(-0.033*(k+2)*h)-21*exp(-0.04*(k+2)*h);
deltaU(k) = deltaU(k-2)*( y(k+2)-2*y(k+1)+y(k) )/(x(1,k)-2*x(1,k-1)+x(1,k-2) );
u(k) = u(k-2) +deltaU(k) ;

u(k+1) = y(k) - x(1,k);
un = u(k);
tt = k*h:h/10 : (k+1)*h;
[T,Y] = ode45(@xslv,tt,x(:,k));
x(:,k+1) = Y(end, : )' ;
end
Y(end) = [];

t=0:h:Tend-1;
t1=0:h:Tend;
t2=0:h:Tend+1;
x1=x(1, : ) ;
x2=x(2, : ) ;
% Выводим результат в виде графиков
figure;
plot(t,deltaU,'k');
xlabel('t, sec'),ylabel('u(t)') ;
legend('u(t)');
grid on;
figure;
plot (t1,x1, 'k', t2,y, 'k--');
xlabel('t, sec'),ylabel ('x(t)') ;
legend('x(t)-real','x(t)-etalon');
grid on;
figure;
plot (t1,x2, 'k', t2,y, 'k--');
xlabel('t, sec'),ylabel ('x(t)') ;
legend('x(t)-real','x(t)-etalon');
grid on;

function dx = xslv(t,x)
global a1 a2 a3 F1 F2 F3 un
% T1=0.5;
% T2=0.25;
% T3=0.05;
% T4=1;
dx = zeros (3,1);
dx(1)=x(2);
dx(2)=x(3);
dx(3)=F3*un-a1*x(3)-a2*x(2)-a3*x(1);

Помогите пожалуйста