Решение систем дифференциальных уравнений методом преобразования Лапласа

Решение систем дифференциальных уравнений методом преобразования Лапласа

Министерство образования и науки Российской Федерации

ФГАОУ ВО "Севастопольский государственный университет"

Институт Информационных технологий и управления в технических системах

Кафедра Информатики и управления в технических системах

Пояснительная записка

к курсовой работе

по дисциплине "Специальные разделы математики"

на тему Решение систем дифференциальных уравнений

методом преобразования Лапласа

Работу выполнил

студент гр. УТС/б-31з К.С. Каратаев

Севастополь - 2015

Содержание

Введение

. Цель курсовой работы

. Задание на курсовую работу

. Решение систем линейных стационарных ОДУ методом преобразования Лапласа

. Постановка задачи и решение

. Стабилизация движения путем введения отрицательной обратной связи

Выводы

Список литературы

Приложение

Введение

Синтез САУ начинается с изучения управляемого объекта и формулирования требований к системе. В соответствии с постановкой задачи из анализа математической модели объекта определяют его программные движения (в частности, состояния равновесия). В реальных условиях программные движения абсолютно точно выполнить невозможно. Поэтому следующим этапом является построение математической модели управляющей системы, обеспечивающей при наличии начальных отклонений и внешних воздействий выполнение программы с необходимой точностью. Синтезируемая модель должна быть устойчивой и удовлетворять требованиям качества переходных процессов. Кроме того, эта модель должна быть физически реализуема с применением элементов, отвечающих требованиям стоимости, надежности, специфическим условиям работы системы и т. п.

При практическом синтезе систем управления широкое распространение получил метод квадратичной оптимизации. Достоинства этого метода заключаются, прежде всего, в линейной структуре регулятора, обеспечивающей простоту анализа промежуточных результатов и реализации системы, а также в простоте вычислительных процедур, опирающихся на развитый пакет программного обеспечения для синтеза линейных квадратичных гауссовских (ЛКГ) регуляторов на ЦВМ.

Линейно-квадратичное гауссовское (ЛКГ) управление относится к современным методам управления. Методология синтеза контроллера позволяет отнести системы управления, построенные на таком принципе, к оптимальным системам, в которых оптимизация проводится по некоторому заданному квадратичному критерию качества. Также эта теория принимает в расчёт присутствие возмущений в виде гауссовабелого шума. Однако несмотря на то, что синтез ЛКГ-контроллеров предусматривает систематическую процедуру расчёта для оптимизации качества системы, главным его недостатком является то, что в рассмотрение не принимается робастность системы. Поэтому ЛКГ-синтез проводится только для систем, имеющих надёжную и точную линейную динамическую модель. Для повышения робастности системы управления применяют более сложные алгоритмы, такие как минимаксный ЛКГ синтез, или комбинированный ЛКГ/H∞синтез. ЛКГ контроллеры могут использоваться как для дискретных, так и для непрерывных систем.

дифференциальный лаплас стабилизация

1. Цель курсовой работы

Целью курсовой работы является овладение умением и практическими навыками решения систем обыкновенных дифференциальных уравнений методом преобразования Лапласа и освоение вычислительных процедур с помощью математических пакетов MathCAD или MatLab в применении к задаче синтеза линейно-квадратичного регулятора с внешними возмущениями в виде гауссового белого шума(LQG).

. Задание на курсовую работу

В качестве объекта исследований выбрана система обыкновенных дифференциальных уравнений (ОДУ), описывающая боковое перемещение самолета относительно заданного курса

В этих уравнениях:

-Y(t) - боковое отклонение: отклонение центра масс самолёта в горизонтальной плоскости относительно заданного курса (в метрах);

-f(t) - угол курса: угол между направлением заданного курса и текущим направлением вектора скорости самолёта (в градусах);

-g(t) - угол крена самолёта (в градусах);

-u(t) - внешнее воздействие: текущее значение угла поворота элеронов самолёта (в градусах), с помощью которого можно управлять перемещением самолёта.

Согласно заданию (вариант №2) коэффициенты модели и начальные условия имеют следующие значения:

= 2,6 = 40

= 0,01 = 0

= 0,05 = -1

= 0,05

В задаче необходимо рассмотреть два вида внешних воздействий u(t) системы ОДУ (2.1) - ступенчатое и гармоническое .

Для данного варианта

= 1(t-30)-1(t-10)

=

Уравнение (2.1) описывает боковое перемещение нестабилизированного самолёта, так что его собственное движение является неустойчивым. Для стабилизации самолёта необходимо сформировать управление в виде отрицательной обратной связи (ООС) по состоянию u(t)=uoc(t), которое может быть представлено в виде

= (2.2)

Коэффициенты обратной связи ky, и kg следует выбрать таким образом, чтобы собственные числа системы уравнений (2.1) и (2.2), описывающей движение стабилизированного самолёта, были равны заданным :

;

В курсовой работе необходимо, с помощью преобразования Лапласа, провести исследование решений систем дифференциальных уравнений для нестабилизированного самолёта (2.1) и для стабилизированного самолёта (2.1), (2.2).

. Решение систем линейных стационарных ОДУ

Методом преобразования Лапласа

Наиболее эффективным методом интегрирования систем линейных ОДУ с постоянными коэффициентами считается метод преобразования Лапласа.

Будем считать, что система ОДУ задана в векторно-матричной форме

Напомним, что здесь x(t) - n-вектор искомых функций, u(t) и f(t) - m- и n-векторы внешних воздействий, соответственно, управляющего и независимого. A(nn) и B(nm) - числовые матрицы коэффициентов системы ОДУ. Предполагается, что для системы (2.3) известны значения искомого вектора x(t) в начальный момент времени t=0, т.е., известны начальные условия

Алгоритм решения (интегрирования) системы линейных стационарных ОДУ с использованием метода преобразования Лапласа содержит три основных шага.

1.Применение преобразования Лапласа к системе (2.3);

2.Решение системы линейных алгебраических уравнение в области изображений и нахождение изображения искомого вектораX(s);

3.Получение оригинала искомого вектора x(t) с помощью обратного преобразования Лапласа.

Выполнение шага 1. Применим преобразование Лапласа к системе ОДУ (2.3) и учтём свойства линейности (1.) и дифференцирования оригинала (2.).

Ls·X(s)-x0=A·X(s)+B·U(s)+F(s) s·X(s) - A·X(s)=B·U(s)+ F(s) + x0 (s·I - A)·X(s) = B∙U(s) + F(s) + x0. (2.4)

Здесь X(s) ÷ x(t), U(s )÷ u(t), F(s ÷ f(t), x0 - вектор начальных значений, I - единичная матрица размерности n.

Выполнение шага 2. Из последнего соотношения получим явное выражение для изображения X(s) искомой вектор-функции x(t)

X(s) = ∙ [B∙U(s)+F(s)] +x0 =(s)+ (s) (2.5)

Первое слагаемое в (2.5)

(s) =∙ [B∙U(s) + F(s)] (2.6)

представляет собой изображение вынужденной составляющей общего решения системы ОБУ.

Второе слагаемое в (2.7)

(s) = xo (2.7)

является изображением собственной составляющей общего решения системы ОДУ.

В силу свойства линейности преобразования Лапласа в изображении вынужденного движения (2.6) можно выделить два слагаемых XB(s)= XU(s)+XF(s) :

-изображение вынужденного движения, вызванного управляющим внешним воздействием u(t)

=; (2.8)

-изображение вынужденного движения, вызванного независимым (неуправляемым) внешним воздействием f(t)

. (2.9)

Выполнение шага 3. Применяя обратное преобразование Лапласа к изображениям XC(s), XU(s) и XF(s), определим функции-оригиналы всех трёх видов движений xC(t), xU(t) и xF(t), сложив которые, получим вектор-функцию x(t)=xC(t)+xU(t)+xF(t) - искомое решение системы линейных ОДУ с постоянными коэффициентами (2.3).

Как следует из формул (2.7), (2.8) и (2.9), изображения XC(s), XU(s) и XF(s) составляющих решения (движения) содержат в качестве первого множителя матричную функцию (s∙I-A)-1 комплексной переменной s.Структурно матричную функцию (s∙I-A)-1 можно представить в известной форме

(2.10)

Здесь det(s∙I-A) - определитель (детерминант) матричной функции (s∙I-A), представляющий собой полином n-ой степени с вещественными коэффициентами ai, i=0,1,... ,n-1 относительно комплексной переменной s - характеристический полином матрицы A

Числовые вещественно-значные матрицы Pn-1, Pn-2,..., P1, P0 служат матричными коэффициентами, образующими матричный (n-1)-степени полином переменной s - присоединённую матрицу adj(s∙I-A).

Нули характеристического полинома (2.11) вместе с нулями знаменателей дробно-рациональных изображений внешних воздействий U(s) и F(s) составляют всё множество полюсов изображения X(s) искомой вектор-функции x(t).

Используя пятое свойство преобразования Лапласа - образуя сумму вычетов функции-изображения X(s)∙est в этих полюсах - получаем окончательное выражение для функции-оригинала x(t) - решение системы ОДУ (2.3).

. Постановка задачи и решение

Дана система из трёх линейных ОДУ с постоянными коэффициентами в нормальной форме (в форме Коши), описывающая управляемое боковое перемещение нестабилизированного самолёта (ЛА - летательного аппарата) относительно заданного курса

(2.12)

В этих уравнениях введены следующие обозначения:

-Y(t) - боковое отклонение: отклонение центра масс самолёта в горизонтальной плоскости относительно заданного курса (в метрах);

-f(t) - отклонение угла курса: угол между направлением заданного курса и текущим направлением вектора скорости самолёта (в градусах);

-g(t) - отклонение угла крена самолёта (в градусах);

-u(t) - внешнее управляющее воздействие: текущее значение угла поворота элеронов ЛА (в градусах), с помощью которого можно управлять боковым перемещением самолёта;

-aYg, afg, agg, agu - коэффициенты математической модели движения ЛА.

Известны значения переменных модели Y(t), f(t), g(t) в начальный момент времени t0=0 - начальные условия - Y0=Y(0), f0=f(0), go=g(0).

Согласно заданию:

= 2,6 = 40

= 0,01 = 0

= 0,05 = -1

= 0,05

Предполагается, что внешнее управляющее воздействие u(t) может быть двух видов - ступенчатое u1(t) и гармоническое u2(t). В качестве ступенчатого воздействия выбирается прямоугольный импульс единичной амплитуды и длительностью 30 секунд - u1(t)=1(t-30)-1(t-10). Здесь 1(t) представляет функцию единичного скачка - функцию Хэвисайда. Гармоническая функция (t)=cos(0,25t) - синусоида единичной амплитуды с частотой 0,25 рад/сек и нулевой фазой.

Выходной переменной (управляемым выходом) определим боковое отклонение Y(t) центра масс ЛА относительно заданного курса полёта. Поведение объекта управления (системы ОДУ) рассматривается на интервале изменения аргумента t (времени) от 0 до 100 секунд.

При подстановке конкретных значений коэффициентов система (2.12) примет вид

(2.13)

Приведем исходную систему (2.12), (2.13) к виду модели в форме пространства состояний

Введём необходимые векторы и матрицы

; ; ; (2.14)

; ; ;

В рассматриваемом примере (2.13), (2.14): искомая вектор-функция x(t) имеет размерность 3; внешнее управляющее воздействие u(t) имеет размерность 1 (скалярная функция); внешнее неконтролируемое воздействие f(t) отсутствует; выходная управляемая переменная Y(t)=x1(t) - скалярная функция; числовая матрица коэффициентов А - квадратная, размерности 3x3; числовая матрица коэффициентов при внешнем воздействии B - вектор столбец размерности 3x1; числовая матрица выхода C - вектор-строка размерности 1x3; управляющее воздействие u(t) непосредственного влияния на выходную переменную Y(t) не оказывает, поэтому матрица D - нулевой скаляр.

Первый шаг алгоритма решения системы (2.13). Применяя преобразование Лапласа к системе ОДУ (2.13) с учётом обозначений (2.14), получим выражение (2.15) для нашего примера (при условии f(t)≡0) в следующем виде

(2.15)

Второй шаг алгоритма. Получив изображения по Лапласу двух видов внешних управляющих воздействий ступенчатого и гармонического (при получении изображения используем третье свойство преобразования Лапласа)

запишем явное выражение для изображения X(s) общего решения в виде

(2.16)

Изображения двух вынужденных движений XS(s) и XG(s) в форме (2.6) для двух видов управляющих воздействий, соответственно ступенчатого и гармонического

,

.

Изображение XC(s) собственного движения в виде (2.7)

.

Третий шаг алгоритма. Находим характеристический полином матрицы коэффициентов A

определяем присоединённую матрицу adj(sI-A)

adj(sI-A)

и формируем матрицу-резольвенту

.

Получим в явном виде изображения собственного и вынужденных управляемых и движений

Оригиналы трех составляющих движений - собственного и двух вынужденных получим, применив процедуру обратного преобразования Лапласа к соответствующим изображениям (s):

(13*t)/25 - (52*exp(t/20))/5 + 252/5

/5 - exp(t/20)/5

exp(t/20)

)

Общее решение системы (2.13) находится как сумма собственного и вынужденного движений . Выполнив суммирование и приведя подобные члены, получим искомые вектор функции в виде

Построим график этих решений. Из графиков следует, что исследуемый объект является неустойчивым, т.к. с увеличением времени моделирования первые компоненты вектор функций решений неограниченно увеличиваются по модулю.

Рис. 1

5. Стабилизация движения путем введения отрицательной обратной связи

Проведём стабилизацию самолёта. Для этого спроектируем стабилизирующий регулятор с помощью введения отрицательной обратной связи по компонентам вектора состояния в виде дополнительного сигнала uoc(t) на вход объекта управления

Уравнение состояния (2.13) с учётом введённой обратной связи изменится следующим образом

После раскрытия скобок и приведения подобных слагаемых в последнем уравнение системы, получим

Матрица коэффициентов А для полученной системы ОДУ примет вид матрицы системы с обратной связью

Объект управления (2.3), (2.4) называют стабилизированным (асимптотически устойчивым) если все собственные числа матрицы А его уравнения состояния (2.3) имеют строго отрицательные вещественные части (расположены строго в левой полуплоскости плоскости комплексных чисел).

По заданию курсовой работы необходимо определить неизвестные коэффициенты kY, kf и kg так, чтобы матрица Аос имела три известных собственных числа (имела заданный спектр). В нашем случае - это три отрицательных числа [-0,02; -0,02; -0,05]. Такой набор чисел соответствует нулям нормированного полинома

(2.18)

С другой стороны, известно, что собственные числа матрицы являются нулями её характеристического полинома или корнями её характеристического уравнения, поэтому для матрицы Аос получим её характеристический полином

(2.19)

Из требования выбора коэффициентов kY, kf и kg такими, чтобы полином имел бы три нуля, равных -0,02 и -0,05 следует, что у полиномов (2.18) и (2.19) должны быть равными коэффициенты при одинаковых степенях переменной s:

,09 = 0,05( 0,024 = 0,005 0,000020 = 0,0013

.0154 4.8000 2.8000

Подставив значения коэффициентов обратной связи в выражение (2.17), получим матрицу динамики стабилизированного самолета

Вычислим резольвенту матрицы (см. программное решение)

Получив для стабилизированного самолета изображения собственного и применив к каждому из этих изображений обратное преобразование Лапласа, определим оригиналы трех составляющих движений - собственного и двух вынужденных и , а затем запишем выражения для общих движений стабилизированного самолета при двух видах управляющего воздействия и (см. программное решение)

Выводы

Построим графики движений и . Из графиков следует, что исследуемый объект является устойчивым. Из графика второго окна видно, при отсутствии воздействия значения компонент стремятся к нулю.

Однако экстремумы показывают, что значения функций f(t) - угловое отклонение от заданного курса и g(t) - угол крена во время воздействия отклоняются от нулевых почти на 20% , что может приводить к значительным флуктуациям курсовой стабильности.

Список литературы

1.В.А. Крамарь, В.А .Карапетьян Методические указания к выполнению курсовой работы по дисциплине «Специальные разделы математики» Севастополь. СевГУ. 2015.-30с.

2.Андык В.С. Теория автоматического управления. Учебное пособие практическим занятиям. Томск. Издательство ТПУ 2004.-108с.

.Бесекерский В.А., Попов Е.П. Теория систем автоматического регулирования. М. Наука, 1972.-762с.

.Дорф Р. Бишоп Р. Современные системы управления |Перевод с английского Б.И. Копылова| М. Лаборатория базовых знаний.

.Линейно-квадратичное гауссовское управление

.Преобразовамние Лапламса

Приложение

Выполнение задания в пакете MatLab

;

%Исходные данные модели ЛА

A = [0 2.6 0; 0 0 0.01; 0 0 0.05];

B = [0 0 0.05];

C = [1 0 0];

D = 0;

x0 = [40;0;-1];

syms t s;

%Оригиналы и изображения входных воздействий

Us_t = heaviside(t-30)-heaviside(t-10);_t = cos(0.25*t);_s = laplace(Us_t); vpa(Us_s,6);_s = laplace(Ug_t); vpa(Ug_s,6);

%Определение резольвенты матрицы А= s*eye(3)-A; vpa(As,6);

Res = As^-1; vpa(Res,6);_pol = vpa(factor(det(As)),2);

%Определение изображений движений= Res*x0;

Xus_s = Res*B*Us_s; vpa(Xus_s,6);_s = Res*B*Ug_s; vpa(Xug_s,6);

%Определение оригиналов движений_t = ilaplace(Xc);

Xus_t = ilaplace(Xus_s); vpa(simplify(Xus_t),6);_t = ilaplace(Xug_s); vpa(simplify(Xug_t),6);

Xcs_t = Xc_t+Xus_t;_t = Xc_t+Xug_t;

%Построение графиков движений нестабилизированного самолета

tk = 100;

time = 0:0.001:tk;(1,3,1);(Xcg_t(1),time); grid on(x_1(t))(1,3,2);(Xcg_t(2),time); grid on(x_2(t))(1,3,3);(Xcg_t(3),time); grid on

title(x_3(t))

%Построение модели стабилизированного самолета

syms Ky Kf Kg

Koc = [Ky Kf Kg];

Aoc = A-B*Koc; vpa(Aoc,4);= collect(det(s*eye(3)-Aoc),s); vpa(pol,8);= poly([-0.02 -0.02 -0.05]); vpa(koef,8);= acker(A,B,[-0.02 -0.02 -0.05]);= A-B*K;(Aoc);_s=s*eye(3)-Aoc; vpa(Aoc_s,6);_oc=Aoc_s^-1; vpa(Res_oc,6);

%Определение изображений движений_oc = Res_os*x0;

Xus_s_oc = Res_oc*B*Us_s; vpa(Xus_s_oc,6);_s_oc = Res_oc*B*Ug_s; vpa(Xug_s_oc,6);

%Определение оригиналов движений_t_oc = ilaplace(Xc_oc);

Xus_t_oc = ilaplace(Xus_s_oc); vpa(Xus_t_oc,6);_t_oc = ilaplace(Xug_s_oc); vpa(Xug_t_oc,6);

Xcs_t_oc = Xc_t_oc+Xus_t_oc; vpa(simplify(Xcs_t_oc),6);_t_oc = Xc_t_oc+Xug_t_oc; vpa(simplify(Xcg_t_oc),6);

%Построение графиков движений стабилизированного ЛА

figure

subplot(1,3,1);(Xcs_t_oc(1),time); grid on(x_1(t))(1,3,2);(Xcs_t_oc(2),time); grid on(x_2(t))(1,3,3);(Xcs_t_oc(3),time); grid on=

2.6000 0

0 0.0100

0 0.0500=

.0500=

0 0= 0=

_t = heaviside(t - 30) - heaviside(t - 10)_t = cos(t/4)= exp(-30.0*s)/s - (1.0*exp(-10.0*s))/s= s/(s^2 + 0.0625)=

[ s, -2.6, 0]

[ 0, s, -0.01]

[ 0, 0, s - 0.05]=

[ 1/s, 2.6/s^2, 0.52/(s^2*(20.0*s - 1.0))]

[ 0, 1/s, 0.2/(s*(20.0*s - 1.0))]

[ 0, 0, 20.0/(20.0*s - 1.0)]_pol = [ 0.05, s, s, 20.0*s - 1.0]=

/s - 13/(25*s^2*(20*s - 1))

/(5*s*(20*s - 1))

/(20*s - 1)=

(0.026*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(s^2*(20.0*s - 1.0))

(0.01*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(s*(20.0*s - 1.0))

(1.0*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(20.0*s - 1.0)=

.01/((20.0*s - 1.0)*(s^2 + 0.0625))/((20.0*s - 1.0)*(s^2 + 0.0625))

Xc_t =

(13*t)/25 - (52*exp(t/20))/5 + 252/5

/5 - exp(t/20)/5

exp(t/20)=

.026*heaviside(1.0*t - 10.0)*(10.0*t - 400.0*exp(0.05*t - 0.5) + 0.5*t^2 + 250.0) + 0.026*heaviside(1.0*t - 30.0)*(10.0*t + 400.0*exp(0.05*t - 1.5) - 0.5*t^2 - 250.0)

.01*heaviside(1.0*t - 30.0)*(20.0*exp(0.05*t - 1.5) - 1.0*t + 10.0) + 0.01*heaviside(1.0*t - 10.0)*(t - 20.0*exp(0.05*t - 0.5) + 10.0)(1.0*t - 30.0)*(exp(0.05*t - 1.5) - 1.0) - 1.0*heaviside(1.0*t - 10.0)*(exp(0.05*t - 0.5) - 1.0)=

.016*cos(0.25*t) - 0.08*sin(0.25*t) + 0.4*exp(0.05*t) - 0.416

.00769231*exp(0.05*t) - 0.00153846*sin(0.25*t) - 0.00769231*cos(0.25*t)

.192308*sin(0.25*t) - 0.0384615*cos(0.25*t) + 0.0384615*exp(0.05*t)_t =

(13*t)/25 - (52*exp(t/20))/5 + (13*heaviside(t - 10)*(20*t - 400*exp(t/20 - 1/2) + (t - 10)^2/2 + 200))/500 - (13*heaviside(t - 30)*(20*t - 400*exp(t/20 - 3/2) + (t - 30)^2/2 - 200))/500 + 252/5

(heaviside(t - 10)*(t - 20*exp(t/20 - 1/2) + 10))/100 - exp(t/20)/5 + (heaviside(t - 30)*(20*exp(t/20 - 3/2) - t + 10))/100 + 1/5

heaviside(t - 30)*(exp(t/20 - 3/2) - 1) - heaviside(t - 10)*(exp(t/20 - 1/2) - 1) - exp(t/20)_t =

(13*t)/25 + (2*cos(t/4))/125 - 10*exp(t/20) - (2*sin(t/4))/25 + 6248/125

/5 - (5*exp(t/20))/26 - sin(t/4)/650 - cos(t/4)/130

(5*sin(t/4))/26 - (25*exp(t/20))/26 - cos(t/4)/26= [ Ky, Kf, Kg]=

[ 0, 2.6, 0]

[ 0, 0, 0.01]

[ -0.05*Ky, -0.05*Kf, 0.05 - 0.05*Kg]

ans = s^3 + (0.05*Kg - 0.05)*s^2 + 0.0005*Kf*s + 0.0013*Ky= [ 1.0, 0.09, 0.0024, 0.00002]= 0.0154 4.8000 2.8000=

2.6000 0

0 0.0100

.0008 -0.2400 -0.0900=

.0200 + 0.0000i

.0200 - 0.0000i

.0500 + 0.0000i=

[ s, -2.6, 0]

[ 0, s, -0.01]

[ 0.000769231, 0.24, s + 0.09]=

[ (20.0*(2500.0*s^2 + 225.0*s + 6.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (1300.0*(100.0*s + 9.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), 1300.0/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)]

[ -0.384615/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (500.0*s*(100.0*s + 9.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (500.0*s)/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)]

[ -(38.4615*s)/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), -(100.0*(120.0*s + 1.0))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0), (50000.0*s^2)/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)]_oc =

(800*(2500*s^2 + 225*s + 6))/(50000*s^3 + 4500*s^2 + 120*s + 1) - 1300/(50000*s^3 + 4500*s^2 + 120*s + 1)

(500*s)/(50000*s^3 + 4500*s^2 + 120*s + 1) - 200/(13*(50000*s^3 + 4500*s^2 + 120*s + 1))

(50000*s^2)/(50000*s^3 + 4500*s^2 + 120*s + 1) - (20000*s)/(13*(50000*s^3 + 4500*s^2 + 120*s + 1))=

(65.0*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)

(25.0*s*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)

(2500.0*s^2*(exp(-10.0*s)/s - (1.0*exp(-30.0*s))/s))/(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0)=

(65.0*s)/((s^2 + 0.0625)*(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0))

(25.0*s^2)/((s^2 + 0.0625)*(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0))

(2500.0*s^3)/((s^2 + 0.0625)*(50000.0*s^3 + 4500.0*s^2 + 120.0*s + 1.0))_t_oc =

(460*exp(-t/50))/9 - (100*exp(-t/20))/9 + (7*t*exp(-t/50))/15

(25*exp(-t/20))/117 - (25*exp(-t/50))/117 - (7*t*exp(-t/50))/1950

(8*exp(-t/50))/117 - (125*exp(-t/20))/117 + (7*t*exp(-t/50))/975=

.0*heaviside(1.0*t - 10.0)*(0.444444*exp(0.5 - 0.05*t) + 0.555556*exp(0.2 - 0.02*t) + 0.0333333*exp(0.2 - 0.02*t)*(t - 10.0) - 1.0) - 65.0*heaviside(1.0*t - 30.0)*(0.444444*exp(1.5 - 0.05*t) + 0.555556*exp(0.6 - 0.02*t) + 0.0333333*exp(0.6 - 0.02*t)*(t - 30.0) - 1.0)

.0*heaviside(1.0*t - 30.0)*(0.0222222*exp(1.5 - 0.05*t) - 0.0222222*exp(0.6 - 0.02*t) + 0.000666667*exp(0.6 - 0.02*t)*(t - 30.0)) - 25.0*heaviside(1.0*t - 10.0)*(0.0222222*exp(0.5 - 0.05*t) - 0.0222222*exp(0.2 - 0.02*t) + 0.000666667*exp(0.2 - 0.02*t)*(t - 10.0))

.0*heaviside(1.0*t - 10.0)*(0.00111111*exp(0.5 - 0.05*t) - 0.00111111*exp(0.2 - 0.02*t) + 0.0000133333*exp(0.2 - 0.02*t)*(t - 10.0)) - 2500.0*heaviside(1.0*t - 30.0)*(0.00111111*exp(1.5 - 0.05*t) - 0.00111111*exp(0.6 - 0.02*t) + 0.0000133333*exp(0.6 - 0.02*t)*(t - 30.0))=

.13944*exp(-0.02*t) - 0.0759527*sin(0.25*t) - 1.11111*exp(-0.05*t) - 0.0283338*cos(0.25*t) - 0.0137785*t*exp(-0.02*t)

.0027244*sin(0.25*t) - 0.00730314*cos(0.25*t) + 0.0213675*exp(-0.05*t) - 0.0140644*exp(-0.02*t) + 0.000105988*t*exp(-0.02*t)

.06811*cos(0.25*t) + 0.182579*sin(0.25*t) - 0.106838*exp(-0.05*t) + 0.0387276*exp(-0.02*t) - 0.000211977*t*exp(-0.02*t)=

.1111*exp(-0.02*t) - 11.1111*exp(-0.05*t) - 65.0*heaviside(1.0*t - 30.0)*(0.444444*exp(1.5 - 0.05*t) - 0.444444*exp(0.6 - 0.02*t) + 0.0333333*t*exp(0.6 - 0.02*t) - 1.0) + 65.0*heaviside(1.0*t - 10.0)*(0.444444*exp(0.5 - 0.05*t) + 0.222222*exp(0.2 - 0.02*t) + 0.0333333*t*exp(0.2 - 0.02*t) - 1.0) + 0.466667*t*exp(-0.02*t)

.213675*exp(-0.05*t) - 0.213675*exp(-0.02*t) - 25.0*heaviside(1.0*t - 10.0)*(0.0222222*exp(0.5 - 0.05*t) - 0.0222222*exp(0.2 - 0.02*t) + 0.000666667*exp(0.2 - 0.02*t)*(t - 10.0)) + 25.0*heaviside(1.0*t - 30.0)*(0.0222222*exp(1.5 - 0.05*t) - 0.0222222*exp(0.6 - 0.02*t) + 0.000666667*exp(0.6 - 0.02*t)*(t - 30.0)) - 0.00358974*t*exp(-0.02*t)

.0683761*exp(-0.02*t) - 1.06838*exp(-0.05*t) + 2500.0*heaviside(1.0*t - 10.0)*(0.00111111*exp(0.5 - 0.05*t) - 0.00111111*exp(0.2 - 0.02*t) + 0.0000133333*exp(0.2 - 0.02*t)*(t - 10.0)) + 0.00717949*t*exp(-0.02*t) - 2500.0*heaviside(1.0*t - 30.0)*(0.00111111*exp(1.5 - 0.05*t) - 0.00111111*exp(0.6 - 0.02*t) + 0.0000133333*exp(0.6 - 0.02*t)*(t - 30.0))=

.2506*exp(-0.02*t) - 0.0759527*sin(0.25*t) - 12.2222*exp(-0.05*t) - 0.0283338*cos(0.25*t) + 0.452888*t*exp(-0.02*t)

.0027244*sin(0.25*t) - 0.00730314*cos(0.25*t) + 0.235043*exp(-0.05*t) - 0.22774*exp(-0.02*t) - 0.00348376*t*exp(-0.02*t)

.06811*cos(0.25*t) + 0.182579*sin(0.25*t) - 1.17521*exp(-0.05*t) + 0.107104*exp(-0.02*t) + 0.00696751*t*exp(-0.02*t)

>>