Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Федеральное агентство по образованию

ТОМСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ СИСТЕМ УПРАВЛЕНИЯ И РАДИОЭЛЕКТРОНИКИ (ТУСУР)

Кафедра промышленной электроники (ПРЭ)

Учебное методическое пособие

ПРОГРАММИРОВАНИЕ

И.М. Егоров

2006

Содержание

• Введение
• Модель цепи в пространстве состояний
• Получение модели цепи в пространстве состояний на основе системы уравнений Кирхгофа
• Пример построения модели цепи в пространстве состояний
• Получение компонентов модели цепи в пространстве состояний на основе матричных операций MathCad
• Аналитическое решение систем линейных дифференциальных уравнений.
• Матричная экспонента
• Некоторые свойства матричной экспоненты
• Матричная экспонента и преобразование подобия
• Собственные числа и собственные вектора матрицы
• Расчет матричной экспоненты на основе преобразования подобия с использованием функций MathCad
• Решение системы дифференциальных уравнений с использованием матричной экспоненты в MathCad
• Собственные числа, колебательный характер переходного процесса и резонансные явления
• Рекомендации по выбору значений параметров элементов схемы
• Расчет реакции схемы на ступенчатое воздействие
• Реакция цепи на одиночный прямоугольный импульс
• Реакция цепи на периодическую последовательность прямоугольных импульсов
• Получение осциллограмм установившегося режима.
• Трассировка графиков
• Задание на курсовое проектирование
• Выбор значений параметров элементов схемы, обеспечивающих колебательный переходный процесс
• Исследование динамики цепи при включении внешнего воздействия воздействии в виде периодической последовательности импульсов единичной амплитуды скважности г = 0.5 с периодом следования Tp = T.
• Исследование установившегося процесса в цепи при воздействии периодической последовательности импульсов
• Приложение А
• Приложение Б
• Приложение В

Введение

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

Одним из основных приемов разработки новых и анализа существующих систем и устройств является их интерактивное моделирование. Средством интерактивного моделирования служит определенным образом построенная вычислительная среда. С помощью этой среды просчитываются различные варианты построения разрабатываемого объекта, выдается результат, как правило, в виде графиков, таблиц, и далее разработчик сам принимает решение о дальнейшем ходе разработки. Такой подход позволяет использовать для синтеза устройств не формализуемые или не достаточно формализуемые критерии качества их функционирования.

В предлагаемом задании на курсовой проект содержится задача с элементами параметрического синтеза простой электрической цепи, решить которую предлагается с использование математического процессора MathCad. При этом необходимо специальным образом построить рабочий лист MathCad с тем, чтобы было легко изменять значения параметров элементов схемы и сразу наблюдать результат такого изменения. Такая технология построения рабочего листа позволяет сформировать в рабочем поле MathCad подобие электронной таблицы. При этом группу операторов задания исходных данных удобно и располагать непосредственно перед операторами вывода конечных результатов расчета с тем, чтобы они одновременно находились в поле зрения оператора - разработчика.

Изменение величин исходных параметров влечет автоматический пересчет всех остальных вычислительных регионов, включая перерисовку имеющихся графиков. Структура рабочего листа MathCad при такой технологии предполагает его разделение на три функциональные области:

· регионы определений функций,

· регионы задания значений параметров

· регионы вызова функций и вывода результатов.

Таким образом, в рабочем пространстве создается программа решения задачи, составленная с использованием специфического языка интерфейса MathCad.

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

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

Практическая часть пособия содержит 100 различных вариантов заданий, из которых следует по указанному правилу выбрать конкретную задачу для курсового проекта. В «Приложении В» приведен реальный пример выполнения и оформления проекта, это вполне работоспособный текст программы, если его без ошибок перенести в MathCad, то программа «оживет» и будет работать.

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

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

Модель цепи в пространстве состояний.

Электрическая цепь, содержащая n реактивных элементов - индуктивностей и / или емкостей - может быть представлена в виде системы n линейных дифференциальных уравнений 1-го порядка с постоянными коэффициентами. Этими дифференциальными уравнениями связываются токи через индуктивности, напряжения на емкостях и возбуждающие сигналы, в роли которых выступают изменяющиеся во времени величины ЭДС источников напряжения и / или величины токов источников тока.

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

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

Получение модели цепи в пространстве состояний на основе системы уравнений Кирхгофа.

Любая математическая модель строится на основе формальной записи физических закономерностей. Для электрической цепи физические закономерности отражаются законами Кирхгофа с использованием т.н. компонентных соотношений.

Напомним, что первый закон Кирхгофа «запрещает» накопление зарядов в узлах схемы: сумма токов втекающих в узел равна сумме токов вытекающих из узла.

Второй закон Кирхгофа представляет собой одну из форм закона сохранения энергии: работа, совершаемая электрическими силами при перемещении заряда по любому замкнутому контуру в цепи, равна работе внешних сил, разделяющих заряды в источниках. Или, что эквивалентно, алгебраическая сумма падений напряжений на участках контура равна сумме алгебраической сумме ЭДС источников, входящих в этот контур.

Компонентные соотношения связывают напряжения на элементе цепи с величиной протекающего через него тока.

Для активных сопротивлений компонентные соотношения отражают известное выражение закона Ома для участка цепи:

, (1.1)

где i - ток через сопротивление;

U - падение напряжения на сопротивлении;

R - величина сопротивления.

Для реактивных элементов (индуктивностей и емкостей) законы связи токов и напряжений носят дифференциальный характер.

Для индуктивности:

(1.2)

где i_L - ток через индуктивность;

U_L - падение напряжения на индуктивности;

L - величина индуктивности.

Для емкости:

U_С - падение напряжения на емкости;

С - величина емкости.

Система уравнений Кирхгофа позволяет связать напряжения и токи в цепи алгебраическими уравнениями. Для линейных цепей, параметры которых R, L и C не зависят ни от величин токов, протекающих через элементы, ни от напряжений на них, эти уравнения имеют линейный вид. Наличие дифференциальных связей (1.2) и (1.3) приводит к тому, что часть алгебраических уравнений превращаются в дифференциальные.

Пример построения модели цепи в пространстве состояний

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Пусть имеется схема цепи с одним внешним воздействием в виде источника ЭДС E(t), показанная на рисунке 1.1 .

В схеме имеется 6 ветвей и 4 узла, после разметки выделено 3 контура.

Из курса ТОЭ известно, что для такой цепи можно составить 6 уравнений: 3 для баланса токов в узлах (I закон Кирхгофа) и 3 для баланса напряжений в контурах(II закон Кирхгофа).

Запишем систему уравнений Кирхгофа для этой цепи с учетом компонентных соотношений.

1) для узла 1:

2) для узла 2:

3) для узла 3: (1.4)

4) для контураI:

5) для контураII:

6) для контураIII:

Система уравнений (4) связывает 9 переменных:

При правильной записи уравнений Кирхгофа получается система независимых уравнений, и можно выразить любые 6 переменных через оставшиеся 3.Поставим задачу выразить переменные через .

Воспользуемся для такого выражения средствами матричной алгебры математического процессора MathCad (см. п.2.4 и п.3.5 в [1]). Представим систему (1.4) в матричной форме:

X•D = 0, (1.5)

где D - прямоугольная матрица 96 (9 столбцов и 6 строк),

X - вектор-столбец из 9 компонент.

Для нашего примера матрица D и вектор X будут иметь вид:

Преобразуем (1.5) к такому виду:

, (1.7)

где матрицы D0 и D1для рассматриваемого примера имеют вид:

В виду того, что исходная система уравнений линейно независима, матрица D0 имеет обратную, поставленная выше задача имеет единственное решение:

, (1.9)

где G = D0^-1•D1 - матрица 6Ч3 (6 строк и 3 столбца)

В нашем случае матрица G имеет вид:

В системе уравнений (1.9) первые два уравнения дифференциальные, остальные - алгебраические. В выражении (1.10) прямоугольниками выделены компоненты матрицы G, представляющие параметры системы дифференциальных уравнений - модель цепи в пространстве состояний. Выпишем в явной форме систему дифференциальных уравнений:

, (1.11)

где

Решение уравнений (1.11) позволяет найти процессы iL(t) и Uc(t). Оставшиеся алгебраические уравнения системы (1.9) дают выражения токов в ветвях схемы в виде линейной комбинации процессов iL(t), Uc(t) и E(t).

Как правило, модель цепи в пространстве состояния помимо дифференциальных уравнений (ДУ), определяющих ее динамику, включает еще и т.н. уравнения наблюдения. Уравнения наблюдения - это выражения токов, протекающих через какие-либо заданные элементы, и/или напряжений между заданными точками схемы.

Пусть в нашем примере наблюдаемыми величинами являются ток через емкость и напряжение на резисторе R4. Согласно компонентным соотношениям

U_R4 = R4•i4 и .

Получим выражения для наблюдаемых переменных U_R4 и i_C, используя выражения для i4 и из (1.9):

(1.12)

или в матричной форме:

(1.13)

Уравнения (1.11 - 1.13) представляют собой полную модель цепи в пространстве состояний. Дальнейшее решение задачи сводится к отысканию решения системы ДУ (1.11).

Получение компонентов модели цепи в пространстве состояний на основе матричных операций MathCad

Для построения модели цепи в пространстве состояний необходимо получить ряд матриц:

A - матрицу системы ДУ;

b - вектор коэффициентов перед функцией внешнего воздействия E(t);

K - матрицу наблюдения, связывающую наблюдаемые переменные с переменными состояния.

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

Первым создается оператор определения матричной функции, задающий матрицу системы уравнений Кирхгофа (1.6):

Далее следуют операторы, определяющие функции формирования подматриц D0 и D1 из (1.8):

В этих операторах использована функция submatrix из библиотеки MathCad. Эти функция выделяет прямоугольный блок заданного размера из матрицы, записанной первым аргументом. В данном случае это матрица, возвращаемая функцией D(L,C,R1,R2,R3,R4).

Затем идет оператор расчета матрицы G из (1.10):

При попытке выполнить этот оператор может возникнуть коллизия из-за вырожденности матрицы D0. Эта ситуация может быть вызвана либо неверной записью системы уравнений Кирхгофа, либо ошибками при формировании исходной матрицы D.

Для проверки правильности выражений, предшествующих получению матрицы D0, можно вывести выражение для ее определителя в символическом виде. Это выражение должно иметь такую структуру, которая не позволяла бы ему обратиться в нуль при какой-либо комбинации значений параметров элементов схемы. В нашем примере выражение для определителя D0 имеет вид:

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

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

За оператором определения функции, вычисляющей матрицу G, следуют операторы задания функций для расчета матрицы A и вектора b модели цепи в пространстве состояний:

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

Завершают последовательность операторов формирования матриц модели операторы определения функций для расчета матрицы наблюдения:

Функции K1 и K2 формируют строки матрицы наблюдения: K1 -

строку для выражения U_R3, K2 - строку для i_С. В этих функциях использован оператор выделения i - го столбца матрицы ^{< i >} . При формировании K1 из матрицы G^T выделяется столбец с номером 5, при формировании K2 из G^T выделяется столбец с номером 1. Напомним, что в MathCad по умолчанию установлена индексация массивов, начинающаяся с нуля.

Двойное применение операции транспонирования матриц связано с тем, что в MathCad есть только оператор выделения столбца матрицы, и нет оператора выделения строки. Поэтому для выделения i-й строки матрицы приходится сначала ее транспонировать, затем выделить в транспонированной матрице i-й столбец и еще раз транспонировать результат.

В функции расчета матрицы наблюдения K применена библиотечная функция stack, осуществляющая объединение однострочных матриц K1 и K2 в единую двухстрочную матрицу K.

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

Неудобство, вызванное длинными списками параметров у вводимых функций, компенсируется тем, что эти функции становятся независимыми друг от друга и впоследствии могут вызываться с набором фактических параметров, определенным непосредственно перед точкой вызова. Функции A(L,C,R1,R2,R3,R4), b(L,C,R1,R2,R3,R4) и K(L,C,R1,R2,R3,R4) будут использоваться в составе вводимых далее функций.

Аналитическое решение систем линейных дифференциальных уравнений

Аналитическое решение дифференциальных уравнений (ДУ) возможно далеко не во всех случаях. Наиболее общие результаты в плане построения аналитического решения получены для систем линейных ДУ.

(2.1)

где

- вектор-функция искомых решений;

, -

- известная вектор-функция внешних воздействий.

Из курса математики известно, что решение неоднородного ДУ (2.1) может быть получено в форме интеграла наложения (или интеграла Дюамеля):

, (2.2)

где Ф(t) - матричная (nЧn) функция времени;

Х0 - вектор начальных условий.

Матричная функция Ф(t) известна в приложениях под многими названиями: переходная матрица системы, импульсная реакция системы, матрица фундаментальных решений системы однородных ДУ и т.п. Системы, представленные математической моделью (2.1), подразделяют в зависимости от поведения Ф(t) на

· устойчивые, если Ф(t) > 0 при t>?;

· неустойчивые, если Ф(t) > ? при t>?;

· физически реализуемые, если Ф(t) ? 0, при t <0;

· физически нереализуемые. если Ф(t) ? 0, при t <0.

Понятие физической реализуемости связано с невозможностью систем реагировать на воздействие до его появления, реальная система не может «предчувствовать» будущие внешние воздействия. Функцию Ф(t) можно представить как «интенсивность» памяти о величине воздействия по прошествии времени t с момента его возникновения. Чем медленнее затухает Ф(t), тем больший вес имеют прошедшие воздействия в текущем Х(t). Таким образом, чтобы исключить инверсию «памяти» в будущее, для физически реализуемых систем требуется дополнительное определение Ф(t) ? 0 для отрицательных t.

В выражении (2.2) первое слагаемое соответствует решению однородного ДУ, т.е. при f(t) ? 0:

при Ф(0) = I, (2.3)

где ,- единичная матрица,

второе - решение неоднородного ДУ

Матричная экспонента

Функция Ф(t) играет важную роль в теории линейных дифференциальных уравнений. Она имеет специальное название - матричная экспонента, и обозначается как e^At. Свое название функция получила из-за ряда аналогий с обычной скалярной экспоненциальной функцией e^at.

Известно, что решением скалярного ДУ первого порядка

при начальном условии y(0) = 1 является экспонента e^at . Сопоставляя этот факт и (2.3), можно отметить формальную аналогию. Эта аналогия распространяется значительно дальше.

Для скалярной экспоненты имеет место разложение в степенной ряд

.

Если в этом ряду заменить скалярную величину а на квадратную матрицу А и единицу на единичную матрицу Е, получим степенной ряд для матричной экспоненты:

(2.4)

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

Некоторые свойства матричной экспоненты

Отметим некоторые свойства матричной экспоненты e^At , вытекающие из ее представления рядом (2.4).

1) Матричная экспонента e^A существует только для квадратных матриц A, т.е. таких матриц, у которых число строк равно числу столбцов.

2) Матричная экспонента e^A - это квадратная матрица тех же размеров, что и A.

3)Дифференцирование матричной экспоненты по t:

, (2.5)

Соотношение (2.3) следует из (2.5) как частный случай при n =1.

4) Интегрирование матричной экспоненты по t:

(2.6)

или (2.7)

Повторное n-кратное интегрирование матричной экспоненты по времени приводит к более общей формуле:

(2.8)

5) Матричная экспонента eA перестановочна с матрицей - аргументом:

Ae^A = e^AA (2.9)

(из курса математики известно, что далеко не всякие матрицы обладают свойством перестановочности)

6) Для матричной экспоненты справедливо соотношение:

e^A•(t1+t2) = e^A•t1 • e^A•t2 (2.10)

7) Если матрица A диагональная, т.е.:

,

то и матричная экспонента тоже диагональная и имеет вид:

. (2.11)

Матричная экспонента и преобразование подобия

Рассмотрим систему линейных уравнений, связывающую два n - мерных вектора x и у:

y = A•x (2.12)

Представим эти уравнения в другой системе координат. Из курса линейной алгебры известно, что переход в иную систему координат аналогичен умножению векторов на матрицу преобразования T. Для взаимнооднозначного преобразования матрица T должна иметь обратную матрицу T^-1. Пусть y1 и x1 - образы векторов y и x в новой системе координат т.е.: y = T•y1 и x = T•x1. Из условий однозначности преобразований следует y1 = T^-1•y и x1 = T^-1•x. Запишем (2.12) через y1 и x1:

T•y1= A•T•x1 или y1= T^-1•A•T•x1 (2.13)

Система уравнений (2.13) отражает те же связи между векторами, что и система (2.12), только в другой системе координат. Матрицы A и Л = T^-1•A•T в математике называют подобными, а само преобразование матрицы T^-1•A•T - преобразованием подобия. Не трудно убедиться, что пара подобных матриц A и Л связана соотношениями:

Л = T^-1•A•T и A = T• Л •T^-1 (2.14)

Преобразование к другой системе координат целесообразно, если в преобразованных координатах система приобретает, в каком-то смысле, более простой вид. Так, например, если удастся найти такое преобразование T, при котором матрица Л оказывается диагональной, то в преобразованных координатах система (2.13) распадется на n независимых линейных уравнений.

Допустим, что найдено такое преобразование T, которое приводит матрицу A к диагональному виду:

. (2.15)

Не трудно убедиться, что это же преобразование приводит к диагональному виду и A^m - любую степень матрицы A:

. (2.16)

Из (2.16) и (2.4) следует, что преобразование T приводит к диагональному виду и матричную экспоненту e^At:

, (2.17)

что позволяет получить выражение для e^At в исходных координатах:

. (2.18)

Выражение (2.18) позволяет вычислить матричную экспоненту в случае, если матрицы T и Л известны. Поиск матриц T и Л сводится к известной в математике задаче о собственных числах и собственных векторах матрицы.

Собственные числа и собственные вектора матрицы

Собственными числами матрицы A, называют числа л, в общем случае комплексные, при которых определитель матрицы

(2.19)

равен нулю.

Иными словами, собственное число л матрицы A должно удовлетворять уравнению

. (2.20)

Если раскрыть определитель, то он превратится в полином относительно л:

, (2.21)

где n - размер матрицы A, а коэффициенты p_k зависят только от элементов матрицы A. Уравнение (2.20) примет иметь вид:

, (2.22)

Таким образом, задача о поиске собственных чисел матрицы размера nЧn сводится к поиску корней полинома степени n.

В общем случае полином (2.21) может быть представлен в виде произведения:

, (2.23)

где л_i - различные корни полинома;

m_i - кратность корня л_i ;

K - число различных корней.

В математике есть т.н. основная теорема алгебры, которая утверждает: всякий полином степени n имеет в поле комплексных чисел ровно n корней, причем каждый корень считается столько раз, какова его кратность. Это означает, что m₁+m₂+…+m_K = n.

Задача поиска корней полинома в аналитическом виде решена лишь для n ? 4, для n > 4 возможен только их численный поиск.

С собственным числом матрицы связано понятие собственный вектор. Собственным вектором матрицы A, соответствующим собственному числу л_i называют вектор t_i, для которого справедливо соотношение:

A•t_i = л_i•t_i , где i = 1…n. (2.24)

Допустим, что матрица A имеет n различных собственных чисел и, соответственно, n собственных векторов. Составим матрицу T, столбцы которой образованы векторами t_i:

и запишем уравнения (2.24) в матричной форме:

(2.25)

Соотношения (2.24) и (2.25) полностью эквивалентны.

Сделаем еще одно предположение: допустим, что собственные вектора линейно независимы. Это допущение справедливо далеко не для всяких матриц, но мы будем считать, что в наших приложениях оно имеет место. Если вектора t_i линейно независимы, то матрица T будет иметь обратную матрицу T^-1. Умножив слева обе части уравнения (2.25) на T^-1 получим:

(2.26)

Таким образом, получаем, что в рамках сделанных нами предположений, матрица A подобна диагональной матрице и преобразование подобия связано с собственными векторами.

На основании вышеизложенного можно предложить один из способов расчета матричной экспоненты e^At по выражению (2.18), используя функции нахождения собственных чисел и собственных векторов из состава встроенных функций математического процессора MathCad.

Расчет матричной экспоненты на основе преобразования подобия с использованием функций MathCad

В составе встроенных функций векторно-матричной алгебры MathCad имеются функции eigenvals(A) и eigenvecs(A).

Функция eigenvals(A), принимая в качестве аргумента квадратную матрицу A, возвращает набор ее собственных чисел, сформированный в виде вектора.

Функция eigenvecs(A), принимая матрицу A, возвращает матрицу собственных векторов T.

Рассмотрим пример использования этих функций. Зададим квадратную матрицу A в виде матричной константы:

Вызовем функцию расчета собственных чисел этой матрицы и поместим результат в вектор Л:

, (2.27)

и выведем значение вектора Л в численной форме:

Вызовем функцию расчета собственных векторов матрицы A и поместим результат в матрицу T:

, (2.28)

и выведем значение матрицы T в численной форме:

Таким образом, встроенные функции MathCad Позволяют получить собственные функции и собственные вектора матрицы A.

Функции eigenvals(A) и eigenvecs(A) используют для расчета собственных чисел и собственных векторов численные методы, что неизбежно вызывает, хотя и не большую, ошибку. Оценим уровень ошибок, вычислив матрицу T^-1•A•T, которая теоретически должна быть диагональной:

Видим, что получившаяся матрица почти диагональная, на ее диагонали, на самом деле, расположены собственные числа матрицы A, Однако, и недиагональные элементы, хотя и малы, но не равны нулю. Это следствие неизбежных ошибок численных методов. Факт возникновения ошибок численных методов обязательно нужно учитывать при возникновении трудно интерпретируемых результатов.

Для вычисления матричной экспоненты e^At на MathCad нам потребуется сформировать диагональную матрицу вида:

, (2.29)

причем с учетом физической реализуемости, нужно доопределить экспоненты в диагонали так, чтобы они были равны нулю при отрицательном времени, т.е. DD(Л, t) должна быть нулевой при t<0.

Сделать это можно, определив функцию пользователя:

EX(a,t) := if(t<0,0,exp(a•t)) (2.30)

Введенная таким образом функция двух аргументов совпадает с e^at только для t ? 0, для t<0 она тождественно равна нулю.

В составе палитры матричных операций MathCad имеется оператор векторизации функций (см.п.2.7 в [1]). В нашем случае мы подвергнем векторизации введенную нами функцию EX(a,t) по параметру a.

Функция - это векторная функция с элементами:

Применив к векторной функции функцию diag из состава встроенных функций MathCad, получим диагональную матрицу вида:

Матрица (2.32) подобна физически реализуемой матричной экспоненте e^At. Остается только воспользоваться обратным преобразованием подобия, чтобы получить выражение для матричной экспоненты в исходной системе координат:

(2.33)

Выражение (2.33) - оператор MathCad, задающий функцию расчета матричной экспоненты (обозначена Ф(Л,t)) , как функцию вектора собственных значений Л и времени t.

Итак, порядок расчета матричной экспоненты можно представить в виде следующей последовательности операций:

Определяем собственные числа матрицы A с помощью функции eigenvals(A).

Определяем матрицу собственных векторов матрицы A с помощью функции eigenvecs(A).

Создаем скалярную функцию - «физически реализуемую» экспоненту EX(a,t) := if(t<0,0,exp(a•t)).

С помощью оператора векторизации превращаем скалярную функцию EX(a,t) в векторную функцию и, используя функцию diag, формируем диагональную матрицу .

Формируем оператор вычисления матричной экспоненты,.:

Решение системы дифференциальных уравнений с использованием матричной экспоненты в MathCad

Рассмотрим следующую последовательность операторов MathCad:

Здесь представлены функции, предназначенные для аналитического решения системы дифференциальных уравнений (1.11) Полагаем, что рассматриваемой последовательности операторов предшествует фрагмент кода, рассмотренный в п.1.3 и функции расчета необходимых матриц введены.

Первым следует оператор определения скалярной «физически реализуемой» экспоненты EX(a,t), за ним идут операторы задания функций расчета вектора собственных чисел Л(L,C,R1,R2,R3,R4) и матрицы собственных векторов T(L,C,R1,R2,R3,R4). Заметим, что у этих функций сохранен тот же список аргументов, что и у функции расчета матрицы системы. Это связано с тем, что в применяемых встроенных функциях eigenvals и eigenvecs в качестве матрицы-аргумента используется результат работы функции A(L,C,R1,R2,R3,R4). У вводимых далее функций список аргументов еще более удлиняется, дополняясь временем t. Последний оператор в этом фрагменте определяет функцию расчета матричной экспоненты Ф(L,C,R1,R2,R3,R4,t).

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

R1= 10 ом; R2 = 400 ом; R3 = 5 ом; R4=800 ом;

L=0.025 гн; C=510^-6 Ф.

При таком сочетании параметров собственные числа матрицы A являются комплексными и образуют комплексно-сопряженную пару:

,

что обусловливает колебательный характер переходного процесса в цепи, наблюдающийся на графиках рис. 2.1.

Функции, представленные на графиках являются решениями однородного ДУ. Составляющие матричной экспоненты Ф(L,C,R1,R2,R3,R4,t)_0,0. и Ф(L,C,R1,R2,R3,R4,t)_1,0. представляют собой колебания iL(t) и Uc(t), соответственно, при начальных условиях iL(0) = 1А и Uc(0) = 0.

Функции Ф(L,C,R1,R2,R3,R4,t)_0,1. и Ф(L,C,R1,R2,R3,R4,t)_1,1.- это те же iL(t) и Uc(t), но при начальных условиях iL(0) = 0 и Uc(0) = 1 V.



Рис. 2.1. Составляющие матричной экспоненты

Собственные числа, колебательный характер переходного процесса и резонансные явления

В случае некратных собственных чисел матрицы A, любая компонента решения системы однородных ДУ x_i(t), представляющая переходный процесс, выражается линейной комбинацией экспонент:

, (2.34)

где л_k - собственные числа матриц A;

q_ik - постоянные числа;

n - порядок системы ДУ, размер матрицы A.

Для всех вариантов схем предлагаемого курсового проекта порядок системы ДУ равен 2.

Собственные числа могут быть вещественными или образовывать комплексно - сопряженные пары.

В случае комплексно - сопряженных собственных чисел л_1,2 переходный процесс образован синусоидами вида:

Ue^Re(л)t sin(Im(л)t +ц),

где л - любое из комплексно - сопряженных чисел л₁ или л₂,

Re(л) и Im(л) реальная и мнимая части л,

U и ц - постоянные величины, не зависящие от л.

Величину щ₀ = Im(л) называют угловой частотой собственных колебаний, величину б = Re(л), при б < 0, называют затуханием переходного процесса. Иногда для характеристики процесса вместо б используют обратную ей величину - постоянную времени ф, имеющую размерность времени:

. (2.35)

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

С круговой частотой щ₀ связана обычная частота f₀ , измеряемая в герцах, и период собственных колебаний T₀, измеряемый в секундах:

(2.36)

Очевидно, что при | Im(л)| = 0 (случай чисто вещественных собственных чисел), колебаний в переходном процессе не будет, и он будет иметь апериодический характер.

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

В дальнейшем нас будут интересовать схемы с резонансом, т.е. схемы с колебательным характером переходного процесса. Этого необходимо достигнуть выбором величин сопротивлений.

Рекомендации по выбору значений параметров элементов схемы

В схемах предлагаемых вариантов заданий содержится одна индуктивность L = 25мГн, одна емкость С = 1мкФ и до пяти активных сопротивлений. Параметры этих элементов следует выбрать таким образом, чтобы цепь обладала резонансными свойствами, иными словами, нужно добиться того, чтобы в цепи имел место колебательный переходный процесс. Для достижения колебательного переходного процесса необходимо:

выбором величин сопротивлений добиться того, чтобы собственные числа матрицы А стали комплексно - сопряженными;

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

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

(2.37)

Величина о численно равна отношению амплитуды переходного процесса в момент времени t+T₀ к амплитуде процесса в момент времени t. Чем больше величина о, тем больше число периодов собственных колебаний будет наблюдаться на графике переходного процесса. При отрицательных Re(л) и |Im(л)| > 0,о> 0, процесс превращается в апериодический и никаких колебаний не наблюдается.

В нашем примере T₀ = 1.04Ч10^-3 с, ф = 0.461Ч10^-3 с, о=0.105. Как видно из графиков, приведенных на рисунке 2.1, такой выбор значений параметров элементов схемы позволяет наблюдать в поле графика 2 периода собственных колебаний.

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

Колебательный характер переходного процесса обусловлен перетеканием энергии от емкости к индуктивности и обратно. Этому перетеканию «мешают» активные сопротивления, которые рассеивают электрическую энергию, необратимо превращая ее в тепло. Выделим на схеме «контур обмена энергией», контур по которому протекает ток между индуктивностью и емкостью. Для схемы рассматриваемого примера такой контур показан на рисунке 2.2.

Для уменьшения потерь энергии необходимо и «расчистить путь» этому току, задав малые значения сопротивлениям, вошедшим в контур. На рис.2.2 это сопротивления R1 и R3.

Величины сопротивлений в других ветвях схемы, которые отходят от «контура обмена» (на рис.2.2 это R2 и R4) , напротив, следует задать большими, чтобы уменьшить их шунтирующее действие.

Этих мер достаточно для того, чтобы цепь обрела колебательный переходный процесс и, следовательно, стала резонансной. Далее следует в интерактивном режиме подобрать соотношения сопротивлений таким образом, чтобы величина о находилась в пределах [0.1; 0.2].

Расчет реакции схемы на ступенчатое воздействие

Решение неоднородного ДУ зависит от характера изменения функции внешнего возбуждения E(t). Рассмотрим частный случай, когда E(t) представляет собой единичную ступенчатую функцию:

Такое изменение E(t) соответствует включению источника ЭДС в момент t = 0, а искомое решение ДУ представляет собой реакцию схемы на это включение. Будем полагать, что до момента t = 0 схема находилась в нулевом состоянии: токи в ветвях схемы отсутствовали, и напряжения на ее элементах были равны нулю.

Воспользовавшись соотношениями (2.2), (2.7) и (2.10), получаем аналитическое выражение для решения неоднородного ДУ:

, (2.34)

где вектор b, согласно (1.11), имеет вид:.

Выражение (2.34) можно представить в иной форме:

, (2.35)

где X_ust - установившееся значение вектора переменных состояния:

, (2.36)

При реализации вычислений по выражениям (2.34 - 2.35) следует учесть, что они верно отражают поведение процесса только для t ? 0, а при t < 0 компоненты получаемого вектора должны быть тождественно нулевыми. Следовательно, использовать единичную матрицу как константу здесь нельзя, требуется ввести дополнительную матричную функцию, которая бы возвращала единичную матрицу при t ? 0 и нулевую матрицу при t < 0. Такую «физически реализуемую» единичную матрицу можно построить, с помощью матричной функции следующего вида:

Здесь задействованы встроенные функции MathCad: функция identity(n), возвращающая единичную матрицу размера n и функция rows(A), возвращающая число строк в матрице - аргументе. Условное выражение if( t < 0, 0, 1 ) возвращает 0 при t < 0 и 1 при t ? 0.

Последовательность операторов MathCad, выполняющих расчет переходного процесса по соотношению (2.35) можно записать в таком виде:

Графики процессов iL(t) и Uc(t) воспроизведены на рисунках (2.2) и (2.3).

Рис. 2.3. Ток через индуктивность при различных значениях R1

Рис. 2.4. Напряжение на емкости при различных значениях R1

Графики на рисунках 2.3 и 2.4 отражают изменение переменных состояния при переходном процессе при различных величинах сопротивления R1: 10 ом, 50 ом и 2 ома. Видно, что с увеличением R1 происходит демпфирование колебаний.

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

Выполним расчет наблюдаемых переменных. Согласно (1.13), связь наблюдаемых величин, в нашем случае это U_R4(t) и i_C(t), с переменными состояния iL(t), Uc(t) и функцией внешнего воздействия E(t) линейна и задается матрицей наблюдения K, функция вычисления которой была введена в п.1.3. Функция вычисления вектора наблюдаемых величин Y на MathCad выглядит следующим образом:

Библиотечная функция stack использована здесь для дополнения вектора состояний X компонентой E(t) - текущим значением функции внешнего возбуждения.

На рисунках 2.5 и 2.6 показаны графики U_R4(t) и i_C(t), соответствующие процессам, приведенным на рисунках 2.3 и 2.4.

Рис. 2.5. Напряжение на резисторе R4 при разных значениях R1

Рис. 2.6. Ток через емкость при разных значениях R1

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

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

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

Реакция цепи на одиночный прямоугольный импульс

Оператор MathCad, определяющий функцию для расчета динамики вектора состояний цепи в ее реакции на одиночный импульс ЭДС длительностью единичной амплитуды выглядит следующим образом:

Здесь использована предварительно определенная функция

Расчет вектора наблюдаемых переменных осуществляется следующей функцией:

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

Рис. 2.7. Динамика вектора состояний при воздействии импульса различной длительности



Рис. 2.8. Динамика вектора наблюдаемых переменных при воздействии импульса различной длительности

Реакция цепи на периодическую последовательность прямоугольных импульсов

Импульсная последовательность характеризуется двумя параметрами - периодом повторения импульсов Tp и скважностью последовательности . Скважность - это безразмерная величина, равная отношению длительности импульсов к периоду их повторения Tp: . Очевидно, что 0 < < 1

Оператор MathCad, определяющий функцию для расчета динамики вектора состояний цепи в ее реакции на импульсную последовательность единичной амплитуды, создаваемую источником ЭДС, выглядит следующим образом:

Здесь Т - период повторения импульсов, - скважность импульсной последовательности.

Приведенный оператор осуществляет расчет решения неоднородного ДУ для момента времени t с нулевыми начальными условиями при t = 0. Этот расчет производится в два этапа. Сначала определяется, сколько целых периодов импульсной последовательности содержится во временном интервале [0; t] и производится вычисление решения на конец последнего полного периода (цикл for в тексте оператора). Затем решение дополняется до искомого времени t.

Оператор определения функции расчета наблюдаемых переменных для этого случая идентичен предыдущим, различие лишь в вызове функции XP, вычисляющей компоненты вектора состояний цепи и функции EIP, рассчитывающей значения ЭДС:

Определенная ранее функция EIP(t,T,) - генератор последовательности прямоугольных импульсов с периодом следования T и скважностью , создается на основе уже введенной функции EI(t, ) - генератора одиночного импульса длительности .

.

где mod(t,T) - встроенная функция MathCad, вычисляющая остаток от деления t на T.

Периодичность EIP по t обусловлена периодичностью величины остатка от деления линейно растущего t на конечную величину T.

На рисунке 2.9 приведены осциллограммы iL(t) и Uc(t) для рассматриваемого примера. На рисунке 2.10 даны осциллограммы переменных наблюдения.

Ток через индуктивность iL(t)

Напряжение на емкости Uc(t)

Рис. 2.9. Компоненты вектора состояния цепи

Напряжение U4(t)



Ток через емкость iC(t)

Рис. 2.10. Компоненты вектора наблюдаемых переменных

На рисунках 2.9 - 2.10 пунктирной линией изображена осциллограмма импульсного сигнала источника ЭДС. Видно, что кривые iL(t) и Uc(t) непрерывны и имеют заметный излом в моменты прихода фронтов импульсов источника ЭДС, ток через емкость в эти моменты имеет мгновенные скачки.

На графиках видно, что после прихода третьего импульса фрагменты осциллограмм начинают практически повторяться, что свидетельствует о наступлении установившегося режима.

Получение осциллограмм установившегося режима

После включения схемы с периодическим источником ЭДС, через некоторое время в цепи возникает установившийся режим. Все токи и напряжения в ветвях приобретают периодический характер, причем для линейных цепей этот период совпадает с периодом повторения импульсов источника ЭДС.

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

Во-первых, компоненты вектора пространства состояний iL(t) и Uc(t) являются решением системы ДУ, и, следовательно, подчиняются соотношению (2.2).

Во-вторых, компоненты вектора состояний непрерывны, никаких скачков и разрывов на их осциллограммах быть не может.

Вследствие периодичности вектора переменных состояния должны совпадать в любые моменты времени, разнесенные на период следования импульсов источника ЭДС. В частности, в установившемся режиме вектора переменных состояния в начале периода X(0) и в конце периода X(Tp) должны совпадать, таким образом, X(0) = X(Tp) = Xp, где Xp - некоторый, пока не определенный, постоянный вектор. Величина Xp - это вектор начальных условий системы неоднородных ДУ, при которых ее решение, начавшись из этого состояния в начале периода, вновь окажется в нем в конце периода.

Для последовательности прямоугольных импульсов с периодом Tp и скважностью вектор Xp можно получить из следующего соотношения:

откуда следует

Можно заметить, что вектор - это значение реакции цепи на одиночный импульс длительности Tp в момент времени Tp. Для расчета такой реакции в MathCad уже была введена функция XI(L,C,R1,R2,R3,R4,t,), поэтому для расчета Xp можно ею воспользоваться, вызвав ее с соответствующими значениями аргументов. Функция вычисления Xp на MathCad будет выглядеть так:

Функция расчета установившихся колебаний компонент вектора состояний цепи будет иметь вид:

Функция mod использована здесь для периодического отображения произвольного значения времени t на интервал [0; T], где T - период повторения импульсов источника ЭДС.

Расчет значений составляющих вектора наблюдения производится следующей функцией:

На рисунке 2.11 показаны временные диаграммы компонент вектора состояния в установившемся режиме для = 0.5 (график представлен сплошной линией), = 0.25(пунктир) и = 0.75 (штриховая линия). Видно, что при = 0.5 колебания имеют максимальную амплитуду.

На рисунке 2.12 приведены графики установившихся колебаний наблюдаемых переменных



Ток через индуктивность

Напряжение на емкости

Отметим, что если компоненты вектора состояний непрерывны, то величина тока через емкость имеет скачки в моменты прихода фронтов импульсов источника ЭДС

Трассировка графиков

Во многих приложениях удобно снимать информацию непосредственно с графиков, измеряя координаты X и Y какой - либо точки. В графической системе MathCad имеются инструментальные средства трассировки построенных графиков.

Кликнув правой кнопкой мыши на поле графика вызываем окно меню, среди пунктов которого есть пункт Trace… Выбрав этот пункт, вызываем окно трассировки (см. рисунок 2.13)



Рис. 2.13. Трассировка графиков

Помещаем курсор мыши в поле графического региона на выбранную линию графика и нажимаем левую кнопку мыши, возникают две штриховые линий, пересекающиеся в указанной точке на кривой. Это своеобразный «прицел», выделяющий точку на линии графика. При этом в полях окна трассировки отображаются координаты точки, в которой находится перекрестье «прицела». При перемещении указателя мыши вблизи кривой перекрестье как бы прилипает к выбранной линии графика, перемещаясь по ее узлам. Фиксируя координаты точек, отображенные в окне, можем определить положение и значение любой точки графика.

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

Задание на курсовое проектирование

Выбор варианта задания

Номером варианта задания служит число, образованное двумя последними цифрами Вашего пароля. Если эти цифры нули, следует выбрать вариант 100.

Будте внимательны при выборе номера варианта! Работы с неправильно выбранным вариантом задания, как правило, не рецензируются и не зачитываются.

1. Построение графического изображения схемы