Математические задачи в электроэнергетике

Математические задачи в электроэнергетике Одним из способов решения задач в энергетики является моделирование. Моделирование – это процесс исследования объектов познания на их моделях; построение и изучение моделей реально существующих предметов, процессов или явлений с целью получения объяснений этих явлений, а также предсказания явлений, интересующих исследователя. Этапы моделирования Первый этап построения модели предполагает наличие некоторых знаний об объекте-оригинале. Модель отображает (воспроизводит, имитирует) какие-либо существенные черты объекта-оригинала. Для одного объекта может быть построено несколько «специализированных» моделей, концентрирующих внимание на определенных сторонах исследуемого объекта или же характеризующих объект с разной степенью детализации. На втором этапе модель выступает как самостоятельный объект исследования. Одной из форм такого исследования является проведение «модельных» экспериментов, при которых сознательно изменяются условия функционирования модели и систематизируются данные о ее «поведении». Конечным результатом этого этапа является множество (совокупность) знаний о модели. На третьем этапе осуществляется перенос знаний с модели на оригинал - формирование множества знаний. Одновременно происходит переход с «языка» модели на «язык» оригинала. Процесс переноса знаний проводится по определенным правилам. Знания о модели должны быть скорректированы с учетом тех свойство объекта-оригинала, которые не нашли отражения или были изменены при построении модели. Четвертый этап - практическая проверка получаемых с помощью моделей знаний и их использование для построения обобщающей теории объекта, его преобразования или управления им. Математическая модель используется в силу следующих причин: Безопасность, низкая стоимость исследования, отсутствие физических ограничений. Математические модели в энергетике Оборудование и системы, используемые для производства и распределения энергии, имеют достаточно полное и точное математическое описание, что позволяет эффективно использовать технологии математического моделирования на различных этапах проектирования и эксплуатации. Уравнения, описывающие любые электромагнитные процессы, известны и неизменны уже почти 150 лет. Это обстоятельство позволяет осуществлять расчет электромагнитных полей в силовых электроустановках с высокой точностью. Программные средства для математического моделирования На сегодняшний день для математического моделирования применяются следующие типы программных средств: Языки программирования общего назначения (C++, Java, Fortran) Специализированные математические пакеты (Mathcad, SPSS, MATLAB) Системы конечно-элементного моделирования (ANSYS, Comsol Multiphysics) Системы псевдофизического моделирования (Electronics Workbench, Multisim, Simulink) MATLAB — сложная современная система, основанная на объектно-ориентированном языке программирования, интегрируемая со средствами псевдофизического моделирования (Simulink) и системами конечно- элементного моделирования. Содержание 1. Создание функции и script-файла 2. Операторы сравнения 3. Логические операторы 4. Cпециальные символы 5. Комплексные числа 6. Векторы 7. Матрицы 8. Графики 9. Интерполяция 10. Решение системы линейных уравнений 11. Система нелинейных уравнений 12. Полиномы 13. Решение линейного уравнения 14. Решение нелинейного уравнения 15. Дифференциал 16. Интеграл 17. Преобразование Лапласа 18. Обработка данных 19. Интегрирование функций. Метод трапеций. 20. Интегрирование. Метод квадратур. 21. If оператор условного перехода 22. For цикл с известным числом повторений 23. While цикл с неизвестным числом повторений 24. Топологические свойства цепей 25. Четырехполюсники 26. Решение дифференциальных уравнений 27. Система дифференциальных уравнений 28. Переходные процессы 29. Моделирование работы двигателя Лабораторные работы 1 лабораторная_работа_1 цепь постоянного тока 2 лабораторная_работа_2 цепь переменного тока 3 лабораторная_работа_3 АЧХ, ФЧХ, переходный процесс 4 лабораторная_работа_4 трансформатор Окна в MATLAB Командное окно (Окно Command Window) Окно Command Window является для пользователя наиболее важным. Посредством этого окна вводятся математические выражения, получаются результаты вычислений, а также выдаются сообщения, посылаемые системой. Математические выражения пишутся в командной строке после знака приглашения ». Каждая команда в командном окне обрабатывается немедленно после нажатия клавиши ENTER, при этом значения выходных параметров выводятся в то же командное окно, а для рисунков открываются графические окна. Рабочая область (Окно Workspace) В процессе работы используются переменные различных типов. Созданные переменные хранятся в специально отведенной области памяти компьютера. Они не исчезают сами по себе, а только при выходе из программы или с помощью специальных команд. При этом переменные (точнее их значения) можно использовать в любом вводимом нами математическом выражении. Окно Workspace (Рабочая область) предоставляет пользователю список всех переменных, хранящихся в рабочем пространстве. Выбрать можно любую переменную, просмотреть ее содержимое или выполнить какие-либо другие действия. Упомянутые действия выполняются посредством контекстного меню (нужно щелкнуть правой кнопкой мыши по имени переменной в списке). Очистить рабочую область можно с помощью меню EDIT / Clear Workspsce. Окно Current Folder (Текущий каталог) Является аналогом известной программы Проводник Окно истории команд (Command History) Все команды, которые набираются в командной строке Command Window (Окно команд), автоматически образуют список, который и выводится в окне Command History (История команд) (Рис. 1). Чем полезен этот список? Если появилась необходимость повторить ранее выполненную команду, ее можно отыскать в списке Command History (История команд) и, дважды щелкнув по ней левой кнопкой мыши, выполнить. Содержимое данного окна не теряется после выхода из системы и выключения компьютера. Очистить список команд можно только с помощью меню EDIT / Clear Command History. MATLAB Все переменные представляются в виде массивов (матриц). MATLAB — сокращение от MATrix LABoratory. Система MATLAB может использоваться в двух режимах: в режиме непосредственного счета (командный режим) и в режиме программирования. Язык системы MATLAB по своей структуре напоминает популярный командный язык Бейсик. С его помощью можно создавать текстовые модули-функции и модули-сценарии (Script-файлы). Файлы, где хранятся такие модули, имеют расширение *.m и называются М-файлами, а находящиеся в них функции – М-функциями. Файл-функция имеет следующую структуру: function var=f_name (Список_параметров_передаваемых_ значений) СОЗДАНИЕ ФУНКЦИИ И SCRIPT-ФАЙЛА к началу Пример: Функция 1: Функция 2: Функция 3: Script-файл: Результат: ЧИСЛА X и х – разные числа = знак присвоения ОПЕРАТОРЫ СРАВНЕНИЯ к началу < меньше <= меньше или равно > больше >= больше или равно == равно ~= неравно ЛОГИЧЕСКИЕ ОПЕРАТОРЫ к началу @@ логическое И || логическое ИЛИ ~ логическое НЕ @ логическое И для массивов | логическое ИЛИ для массивов CПЕЦИАЛЬНЫЕ СИМВОЛЫ к началу : Определение векторов, индексация массивов, определение циклов с известным числом повторений [ ] Определение массивов, определение функций, возвращающих несколько значений ( ) Скобки в математических выражениях, передача аргументов функций … Перенос части команды в следующую строку текста , Разделение строк массива, разделение параметров функции @ Указатель функции a^x Степенная функция x^a Показательная функция sqrt(x) Квадратный корень exp(x) Экспонента log(x) Натуральный логарифм log10(x) Десятичный логарифм abs(x) Модуль fix(x) Отбрасывание дробной части числа floor(x) Округлениие до меньшего целого ceil(x) Округлениие до большего целого round(x) Обычное округление rem(x,y) Остаток от деленияx x на y без учёта знака mod(x,y) Остаток от деленияx x на y с учёта знака sign(x) Знак числа factor(x) Разложение числа x на простые множители sin(x) Синус sinh(x) Синус гиперболический asin(x) Арксинус asinh(x) Арксинус гиперболический cos(x) Косинус cosh(x) Косинус гиперболический acos(x) Арккосинус acosh(x) Арккосинус гиперболический tan(x) Тангенс tanh(x) Тангенс гиперболический atan(x) Арктангенс inf - бесконечность ВВОД КОМПЛЕКСНЫХ ЧИСЕЛ к началу Мнимая единица – i или j >>х=1+i*2 х=1.0000е+000+2.000е+000i ДЕЙСТВИЯ С КОМПЛЕКСНЫМИ ЧИСЛАМИ x*y x-y x^y abs(x) - модуль комплексного числа real(x) – действительная часть imag(x) – мнимая часть angle(x) – значение аргумента комплексного числа conj(x) – комплексно сопряженное число ВЕКТОРА к началу Ввод вектор-строка Вектор-столбец Вектор, состоящий из двух векторов Вектор-строка Вектор-столбец Сокращенное введение вектора, значение элементов которого составляет арифметическую прогрессию СЛОЖЕНИЕ (вычитание) ВЕКТОРОВ ТРАНСПОНИРОВАННИЕ УМНОЖЕНИЕ ВЕКТОРА НА ЧИСЛО УМНОЖЕНИЕ ВЕКТОРОВ С помощью вектора МОЖНО ЗАДАТЬ значения х, для которых можно будет посчитать значения функции y Или Из каждого элемента вычитается 2 Поэлементное умножение векторов Поэлементное деление векторов Поэлементное возведение в степень Требуется вычислить a^(2-x)*sin(x) для значений a=2, x=0,1,2,3,4,5 МАТРИЦЫ к началу ВВОД Обращение а(2,3) дает значение 6 (2 строка, 3 столбец) Матрица с нулевыми элементами Можно сначала задать матрицу с нулевыми элементами, а потом её заполнить значениями: (при этом сохраняется последнее введенное значение) ones (3,5) – получим матрицу c единичными элементами eye(3,5) – единицы на главной диагонали rand(3,5) –случайные числа Введем матрицы a и b Синус каждого члена матрицы а Поэлементное произведение Поэлементное деление (члены а делятся на b) Поэлементное деление (члены b делятся на а) Поэлементное возведение в степень Прибавление к матрице числа (вычитание числа, умножение, деление на число аналогично) Возведение в степень Сложение (вычитание) матриц Транспонирование Умножение матрицы на матрицу (количество столбцов 1 матрицы = количеству строк 2 матрицы) Обращение матрицы Деление слева направо / ( 4/2 = 2 ) Справа налево \ ( 4\2 = 0,5 ) Вычисление определителя (метод Гаусса ) Определение характеристического полинома ГРАФИКИ к началу Два графика в одних осях Два графика и более (например, один под другим) Функция задана параметрически Задана функция z(t) = t2 + M(t), где M(t) = 2sin(10t) Надо построить график z(t). Задаем функцию М(t) function Mso=M1(t) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here t=0:0.01:10 Mso=2+sin(10*t) end Строим график z(t) t=0:0.01:10; z=t.^2+M1(t); %функция задана в отдельном файле plot(t,z);axis([0 10 0 20]) Или можно построить график в одной программе: t=0:0.01:10; Mso=2+sin(10*t) z=t.^2+M1(t); plot(t,z);axis([0 10 0 20]) ВЕКТОРА КОМПЛЕКСНЫХ ЧИСЕЛ В полярных координатах (hold on после построения вектора а) Вектор l задает масштаб по оси ИНТЕРПОЛЯЦИЯ к началу Интерполяция является процессом вычисления (оценки) промежуточных значений функций, которые находятся между известными или заданными точками.   Обзор функций интерполяции Функции Описание griddata Двумерная интерполяция на неравномерной сетке. griddata3 Трехмерная интерполяция на неравномерной сетке. griddatan Многомерная интерполяция (n >= 3). interp1 Одномерная табличная интерполяция. interp2 Двухмерная табличная интерполяция. interp3 Трехмерная табличная интерполяция. interpft Одномерная интерполяция с использованием быстрого преобразования Фурье. interpn Многомерная табличная интерполяция. pchip Кубическая интерполяция при помощи полинома Эрмита. spline Интерполяция кубическим сплайном. Пример Есть некоторая функция z = t2 + M(t), где t – время, М(t) – величина, которая тоже зависит от времени. Зависимость М от t: t 0.1 0.2 0.3 0.4 0.5 1 1.5 2 3 4 М 1 1 6 6 1 1 1 1 1 1 1 function Mso=M(t) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here x=[0 0.1 0.2 0.3 0.4 0.5 1 1.5 2 3 4]; y=[1 1 6 6 1 1 1 1 1 1 1]; Mso=interp1(x,y,t) end t=0:0.01:10 z=t.^2+M(t) plot(t,z);axis([0 3 0 10]) РЕШЕНИЕ СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ к началу Система линейных алгебраических уравнений (СЛАУ) может быть записана в матричной форме: матрица системы (коэфиициенты СЛАУ) Вектор-столбец неизвестных Вектор-столбец свободных членов А*Х = В Основная функция для решения СЛАУ Х=linsolve(A,B) СЛАУ можно также решить, используя деление Х=А\В СИСТЕМА НЕЛИНЕЙНЫХ УРАВНЕНИЙ к началу Пример №1: Функция solve()  находит аналитическое решение и для хранения результата создаётся структура с полями, соответствующими именам неизвестных. В этом примере имя структуры s , в поле x хранятся значения неизвестной x, в поле y − значения неизвестной y. Эти значения можно вывести: Пример №2: s=solve('x=y^2+2*x-2','sin(x)-y=0') s = x: [1x1 sym] y: [1x1 sym] >> s.x ans = 1.1596952618411803821679304065971 >> s.y ans = 0.9166813722110969228834223995511 Программа: clc s=solve('x=y^2+2*x-2','sin(x)-y=0') s.x s.y ПОЛИНОМЫ к началу p(x)=anxn + … + a2x2 + a1x + a0 Ввод a4 = 1, a3 = 2, a2 = 3, a1 = 4, a0 = 5 p(x)=1x4 + 2x3 + 3x2 + 4x + 5 Умножение полинома р1 на р2 Деление полинома р1 на р2 Найти корни полинома p(x)=x5 + 8x4 + 31x3 + 80x2 +94x + 5 Построение полинома, если известны его корни РЕШЕНИЕ ЛИНЕЙНОГО УРАВНЕНИЯ к началу 2y – 2 = 0 Полином р = 2y – 2 РЕШЕНИЕ НЕЛИНЕЙНОГО УРАВНЕНИЯ Содержание 4у2 – 4у + 1 = 0 Полином р = 4у2 – 4y + 1 или Дифференциал к началу Найти (sin(х))’ Интеграл к началу Найти Преобразование Лапласа к началу ОБРАБОТКА ДАННЫХ к началу Получена зависимость у(х) у 1 3 5 6 8 х 41 2 8 4 5 Зависимость можно задать как матрицу xydata, содержащую 2 строки х и у Определение количество строк, столбцов Интегрирование функций. Метод трапеций. к началу Проинтегрировать функцию Интегрирование. Метод квадратур. к началу Проинтегрировать функцию Пример: ПРОГРАММИРОВАНИЕ if ОПЕРАТОР УСЛОВНОГО ПЕРЕХОДА к началу Алгоритм Синтаксис: If условие 1 операторы 1 elseif (этих блоков может быть несколько) условие 2 операторы 2 else операторы 3 end Пример clear all % очищаем память clc % очищаем экран a=input('загадайте число a = ')% диалоговое окно x=randi([1,2]) % генерация случайного числа из диапазона (1,2) if a~=x disp('не угадал'); else disp('угадал') end Моделирование цепи RLC % задаются параметры цепи R=100; L=0.1; C=0.001; f=50; % расчет xL=2*pi*f*L xC=1/(2*pi*f*C) phaza=atan((xL-xC)/R) t=0:0.000001:0.02; % задается напряжение цепи u=10*sqrt(2)*sin(2*pi*f*t); % вычисляется ток цепи i=1000*sqrt(2)/sqrt(R^2+(xL-xC)^2)*sin(2*pi*f*t+(0-phaza)); % строится график напряжения plot(t,u,'g-');grid hold on % строится график тока plot(t,i,'r-');grid % проверяется условие наличия резонанса if abs(xL-xC)<0.5 disp('резонанс'); elseif abs(xL-xC)<10 disp('почти резонанс'); else disp('резонанса нет'); end Лабораторная работа № 2 к началу Моделирование однофазной цепи переменного тока Вводятся следующие параметры цепи: Программа должна указать, есть ли в цепи резонанс (условие наличия резонанса: если (xL – xC) < 0.5), резонанс есть). Также программа должна рассчитать сопротивление катушки, сопротивление конденсатора, полное сопротивление цепи, действующее значение тока в цепи. Построить зависимости мгновенных значений напряжения u и тока i от времени. Построить векторную диаграмму: вектор тока и векторы напряжения: входное напряжение в цепи Uвх, напряжение на катушке UL, напряжение на конденсаторе Uc, напряжение на резисторе UR. clc % задаются параметры цепи U=220; R=30; L=0.05; C=0.00002; f=50; % расчет параметров xL=2*pi*f*L %считается сопротивление катушки xC=1/(2*pi*f*C)%считается сопротивление конденсатора I=U/(R^2+(xL-xC)^2)%считается действующее значение тока phaza=atan((xL-xC)/R)*180/pi %угод сдвига фаз считается в градусах % проверяется условие наличия резонанса if abs(xL-xC)<0.5 disp('резонанс'); else disp('резонанса нет'); end %теперь строим графики - синусоидальные и вектора % 1. строим синусоидальные графики тока и входного напряжения t=0:0.000001:0.04; % задаем интервал времени от 0 до 0.04 с шагом 0.000001 % задается напряжение цепи u=U*sqrt(2)*sin(2*pi*f*t); % вычисляется ток цепи i=U*sqrt(2)/sqrt(R^2+(xL-xC)^2)*sin(2*pi*f*t+(0-atan((xL-xC)/R))); % строится график напряжения subplot(2,1,1); % будем строитьв одном окне 2 графика. в верхнем окне г график синусоидального напряжения и тока plot(t,u,'g-'); gtext('график напряжения') hold on % команда строит поверх уже построенного графика напряжения график тока % строится график тока plot(t,20*i,'r-');grid; %домножаем ток на 20, чтобы на графике его было видно..т.к. он получается очень маленький axis([0 0.05 -450 450])%задаются значения по осям "х" и "y" gtext('график тока') % построения векторов subplot(2,1,2);%во втором окне строим векторную диаграмму % но сначала надо рассчитать комплексные значения тока и напряжений для j=sqrt(-1)%вводим мнимую единицу j, т.к. i - это мгновенное значение тока I=U/(R+j*xL-j*xC) UR=I*R UL=I*j*xL UC=I*(-j*xC) UR+UL+UC %можно сделать проверку 2 закона кирхгофа. если значение получилось "220", значит, все верно считается % строим вектора тока, напряжений hold on compass(U,'r-') grid; gtext('U') compass(UR,'b-') gtext('UR') compass(UC,'g-') gtext('UC') compass(UL,'k') gtext('UL') compass(10*I,'y-.') gtext('I') for ЦИКЛ С ИЗВЕСТНЫМ ЧИСЛОМ ПОВТОРЕНИЙ к началу Синтаксис for переменная=выражение операторы end Пример 1: программа считает для каждого а значение b for a=[1,2,3] b=a+1 end Пример 2: Исследование режимов работы цепи постоянного тока Определить показание амперметра № U R1 R2 Сопротивление ветви с резистором R3 I 1 10 5 10 Ключ К1замкнут, К2 разомкнут 10 2 10 5 10 Ключи К1и К2 замкнуты 3 10 5 10 Ключ К1 зразомкнут Программа скрин: текст parametrs=[10 5 10 10; 10 5 10 0; 10 5 10 inf] for i=1:3 I=parametrs(i,1)/(parametrs(i,2)+parametrs(i,3)/(parametrs(i,3)/parametrs(i,4)+1)); parametrs(i,5)=I end расчет while ЦИКЛ С НЕИЗВЕСТНЫМ ЧИСЛОМ ПОВТОРЕНИЙ к началу Синтаксис while end Пример a=1 Пример 1: Решить уравнение (метод половинного деления) Уравнение: 2х – 6 = 0 1. Задаем функцию 2. Функцию сохраняем в файле equ.m (название функции должно совпадать с именем файла) 3. Зададим диапазон, в котором, как предполагается, находится корень уравнения: х1… х2. -100…100 4. Зададим точность error = 0,001. 5. В цикле будем делить отрезок [х1…х2], потом проверять, в какой половине находится корень уравнения (знаку). Цикл будет повторяться до тех пор, пока модуль значения функции в точке не будет меньше заданной погрешности error. 6. Если задан интервал, в котором корня нет, используется оператор выхода из цикла break (чтобы не было зацикливания). Пример 2 Найти сопротивление конденсатора для ряда частот f=0, 50, 100, 150 Гц Расчет Топологические свойства цепей к началу R1 R2 R3 R4 R5 R6 Е1 Е2 Е3 Матрица соединений A 1 – узел – начало ветви -1 – узел – конец ветви 0 – ветвь не присоединена к узлу ветви у з л ы 1 2 3 4 5 6 1 1 1 1 2 -1 1 -1 3 1 -1 -1 А = G – диагональная матрица проводимостей. G = 1/R. G = E = Система уравнений в матричной форме (А*G*АT)* UUZ = (-A*(G*E+J)) Токи ветвей I=G*(A’*U+E) + J Пример: R1 R2 R3 R4 R5 R6 Е1 Е2 Е3 5 1 10 15 20 25 10 5 15 Матрица соединений 1 – узел – начало ветви -1 – узел – конец ветви 0 – ветвь не присоединена к узлу ветви у з л ы 1 2 3 4 5 6 1 1 1 1 2 -1 1 -1 3 1 -1 -1 ТРЕХФАЗНАЯ СИСТЕМА Матрица соединений A 1 – узел – начало ветви (ток выходит из узла) -1 – узел – конец ветви (ток входит в узел) 0 – ветвь не присоединена к узлу ветви 1 2 3 4 1 узел 1 1 1 -1 Пример № EA EB EC ZL1 ZL2 ZL3 Za Zb Zc ZN 1 220 220e-j120 220e+j120 5 5 5 100+100i 100+100i 100+100i 5 Программа % Задаем значения ЭДС, сопротивлений линий и фаз ea=220+0i; eb=-110-190.526i; ec=-110+190.526i; z0=5+0i; zl1=5+0i; zl2=5+0i; zl3=5+0i; zna=100+100i; znb=100+100i; znc=100+100i; % Матрица соединений A=[1 1 1 -1]; % проводимости ветвей G1=1/(zl1+zna); G2=1/(zl2+znb); G3=1/(zl3+znc); G4=1/z0; % матрица проводимостей G=[G1 0 0 0 0 G2 0 0 0 0 G3 0 0 0 0 G4]; % вектор-столбец ЭДС и источников тока E=[ea;eb;ec;0]; J=[0;0;0;0]; % вычисления Guz=A*G*A'; Jz=-A*(G*E+J); U=Guz\Jz; I=G*(A'*U+E)+J; ufa=I(1,1)*zna; ufb=I(2,1)*znb; ufc=I(3,1)*znc; u0=I(4,1)*z0; % считаем токи ветвей ia=I(1,1) ib=I(2,1) ic=I(3,1) i0=I(4,1) % строим графики hold on o=[ea,eb,ec];compass(o,'k') gtext('UA') gtext('UB') gtext('UC') compass(100*ia,'r-') gtext('ia') compass(100*ib,'b-') gtext('ib') compass(100*ic,'g-') gtext('ic') compass(100*i0,'k') gtext('in')%ток очень маленький, поэтому его не видно практически на векторной диаграмме compass(ufa,'r-.') gtext('ufa') compass(ufb,'b-.') gtext('ufb') compass(ufc,'g-.') gtext('ufc') compass(u0,'k-.') gtext('uNn')%напряжение очень маленькое, поэтому его не видно практически на векторной диаграмме ЧЕТЫРЕХПОЛЮСНИКИ к началу Построение АЧХ, ФЧХ по передаточной функции Пример №1 a=[1]; % вектор коэффициентов числителя b=[0.001 1]; % вектор коэффициентов знаменателя freqs(a,b) Пример №2 Или w=0:30000; r=10; xL=w*0.025i; xC=-i*1./(w*0.0000005); l=(((1./(r+(r.*(r+xL+xC)./(r+r+xL+xC)))).*(r.*(r+xL+xC)./(r+r+xL+xC)))./(r+xL+xC)).*r; subplot(2,1,1) plot(w,abs(l)) subplot(2,1,2) plot(w,atan(imag(l)./real(l))) Еще пример: 1. w=1:1000; r=10 xL=w.*1i; xC=-i*1./(w*0.0001); l=((xC.*r./(xC+r)))./(xL+(xC.*r./(xC+r))); subplot(2,1,1) plot(w,abs(l)) subplot(2,1,2) plot(w,angle(l)) 2. w=1:1000; r=10 xL=w.*1i; xC=-i*1./(w*0.0001); l=((xC.*r./(xC+r)))./(xL+(xC.*r./(xC+r))); subplot(2,1,1) plot(w,abs(l)) subplot(2,1,2) plot(w,atan(imag(l)./real(l))) 3. r=10 L=1; C=0.0001; a=[1] b=[L*C*r L/r 1]; freqs(a,b) Определение полюсов, нулей передаточной функции, коэффициента усиления системы a=[1]; % числитель b=[0.001 1]; % знаменатель [z,p,k]=tf2zpk(a,b) ОПРЕДЕЛЕНИЕ ВЫХОДНОГО СИГНАЛА ПО ПЕРЕДАТОЧНОЙ ХАРАКТЕРИСТИКЕ Пример №1 Комплексная передаточная функция: i – мнимая единица. Пусть входной сигнал Тогда выходной сигнал: Код программы: w=1000 %частота сигнала w H=1/(0.001*w*i+1) % передаточная функция 141*z % Uвх*H(iw) Решение: РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ к началу или у’ производная 1 порядка Dy или у’’ производная 2 порядка D2y или у’’’ производная 3 порядка D3y Пример 1: , начальные условия не заданы , начальные условия заданы: у(0)=0 Пример 2 , у(0)=2 Пример 3 , у(0)=-0,7 Уравнение II порядка Пример 1: y'’ = -a2y, начальные условия: у(0)=1, y’ = 0. Пример 2: %переходный процесс 2 порядка %ввод диф.уравнения в рабочее поле (апериодический процесс, R=100) y=dsolve('D2y*0.000001+Dy*0.0001*100+y=10','y(0)=0','Dy(0)=0') y = 10 - (2*15^(1/2)*sin(250*15^(1/2)*t))/(3*exp(250*t)) - (10*cos(250*15^(1/2)*t))/exp(250*t) ezplot(y,[0 0.1]);axis([0 0.1 0 20])%построение графика апериодический закон изменения %ввод диф.уравнения в рабочее поле (критический режим, R=20) y=dsolve('D2y*0.000001+Dy*0.0001*20+y=10','y(0)=0','Dy(0)=0') y = 10 - (10000*t)/exp(1000*t) - 10/exp(1000*t) ezplot(y,[0 0.02]);axis([0 0.02 0 20]) %ввод диф.уравнения в рабочее поле (колебательный режим, R=2) y=dsolve('D2y*0.000001+Dy*0.0001*2+y=10','y(0)=0','Dy(0)=0') y = 10 - (10*11^(1/2)*sin(300*11^(1/2)*t))/(33*exp(100*t)) - (10*cos(300*11^(1/2)*t))/exp(100*t) ezplot(y,[0 0.02]);axis([0 0.02 0 20]) %для определения функции тока можно взять производную от «у». СИСТЕМА ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ к началу Для решения дифференциальных уравнений вида: , где - вектор-столбец, F(t) = (t, x1, x2, …, xn) – функция-столбец, зависящая от времени и компонент вектора х предусмотрена функция ode45. ПЕРЕХОДНЫЕ ПРОЦЕССЫ к началу Схема 1: Диф. уравнение цепи: 1 способ: u=dsolve('Du*10*0.0001+u=10','u(0)=0') ezplot(u,[0 0.05]); axis([0 0.02 0 10]) если источник питания синусоидальный u=dsolve('Du*10*0.0001+u=10*sin(3000*t)','u(0)=0') ezplot(u,[0 0.05]); axis([0 0.02 -10 10]) 2 способ: Задаем функцию function dUc=ProcessRC(t,Uc) %Параметры цепи R=10; C=0.0001; U=10; dUc=(U/(R*C)-Uc/(R*C)); % дифференциальное уравнение цепи end Программа [t,Uc]=ode45(@ProcessRC,[0 0.005],0); plot(t,Uc) % % %начальное значение напряжения на конденсаторе % %диапазон изменения времени %название функции Если источник синусоидальный function Duc=Process(t, uc) R=10; U=10*sin(3000*t); C=0.0001; Duc=-uc/(R*C)+U/(R*C); end [t,uc]=ode45(@Process,[0 0.02],0); plot(t,uc) Схема 2: Диф. уравнение цепи: 1 способ: u=dsolve('Du*10*0.0001+u=0','u(0)=10') ezplot(u,[0 0.05]); axis([0 0.02 0 10]) 2 способ: Задаем функцию: function dUc=ProcessRC(t,Uc) %Параметры цепи R=10; C=0.0001; dUc=-Uc/(R*C); % дифференциальное уравнение цепи end Схема 3: 1 способ: iL=dsolve('DiL*0.01+iL*10=10','iL(0)=0') ezplot(iL,[0 0.05]); axis([0 0.02 0 3]) 2 способ: function diL=ProcessRL(t,iL) %Параметры цепи R=10; U=10; L=0.01; diL=-iL*R/L+U/L; % дифференциальное уравнение цепи end [t,iL]=ode45(@ProcessRL,[0 0.005],0); plot(t,iL) % % %начальное значение напряжения на конденсаторе % %диапазон изменения времени %название функции Схема 4: В этой схеме или %ввод диф.уравнения в рабочее поле (периодический процесс, R=5) uc=dsolve('D2uc*0.000001+Duc*0.0001*5+uc=10','uc(0)=0','Duc(0)=0') ezplot(y,[0 0.05]);axis([0 0.05 0 20]) Решение и график: uc = 10 - (2*15^(1/2)*sin(250*15^(1/2)*t))/(3*exp(250*t)) - (10*cos(250*15^(1/2)*t))/exp(250*t) Ток iL можно найти, продифференцировав uC и также построить график. il=0.0001*diff(uc) ezplot(il,[0 0.05]);axis([0 0.05 -1 1]) il = (4*15^(1/2)*sin(250*15^(1/2)*t))/(15*exp(250*t)) В этой схеме Или в одной программе: uc=dsolve('D2uc*0.000001+Duc*0.0001*5+uc=10','uc(0)=0','Duc(0)=0') il=0.0001*diff(uc) subplot(2,1,1) ezplot(uc,[0 0.05]);axis([0 0.05 -10 19]),grid subplot(2,1,2) ezplot(il,[0 0.05]),axis([0 0.025 -0.5 0.9]),grid В форме КОШИ Составляем дифференциальные уравнения для переменных состояния uc и iL. Система уравнений в форме Коши В матричном виде: = ∙ + f(1)- это iL, f(2)- это uC. Начальные условия Для решения дифференциальных уравнений в форме Коши предусмотрена функция ode45 текст function dUdt=ProcessRLC(t,f) R=5; L=0.01; C=0.0001; U=10; dUdt=[-R/L*f(1)-1/L*f(2)+U/L;1/C*f(1)]; end Или можно функцию dUdt представить в матричном виде: [t,f]=ode45('ProcessRLC',[0 0.05],[0 0]); subplot(2,1,1) plot(t,f(:,1),'r-'),axis([0 0.025 -0.5 0.9]),grid subplot(2,1,2) plot(t,f(:,2)),axis([0 0.05 -10 19]),grid Если синусоидальный источник function dUdt=ProcessRLC(t,f) R=5; L=0.01; C=0.0001; U=10*sin(3000*t); dUdt=[-R/L*f(1)-1/L*f(2)+U/L;1/C*f(1)]; end [t,f]=ode45('ProcessRLC',[0 0.05],[0 0]); subplot(2,1,1) plot(t,f(:,1),'r-'),axis([0 0.025 -0.5 0.9]),grid subplot(2,1,2) plot(t,f(:,2)),axis([0 0.05 -10 19]),grid НЕЛИНЕЙНАЯ ЦЕПЬ. ПЕРЕХОДНЫЙ ПРОЦЕСС function dUdt=ProcessRLC(t,f) R=5; L=0.01; C=0.0001; U=10; dUdt=[-(R*f(1)^2+1*f(1)^5+10)/L*f(1)-1/L*f(2)+U/L;1/C*f(1)]; end [t,f]=ode45('ProcessRLC',[0 0.05],[0 0]); subplot(2,1,1) plot(t,f(:,1),'r-'),axis([0 0.025 -0.5 0.9]),grid subplot(2,1,2) plot(t,f(:,2)),axis([0 0.05 -10 19]),grid Задача 2 U= 120sin(1178*t). rвх = 100 Ом, rн = 200 Ом, r3 = 10 Ом, C = 2мкФ, L = 0.15 + 0.05m, m = 8, L = 0,55 Гн, r5 = 20 Ом. Вывод формулы На основании первого закона Кирхгофа: i1 = ic + iL+i4; на основании второго закона Кирхгофа: uBx = i1·rBx + uL; -uc - ic·r3 + uL = 0; -uL + i4(r5 + rн) = 0. Метод переменных состояния основывается на двух уравнениях, записываемых в матричной форме. В качестве переменных состояния в электрических цепях следует выбрать токи  в индуктивностях и напряжения  на емкостях, причем не во всех индуктивностях и не на всех емкостях, а только для независимых, таких, которые определяют общий порядок системы дифференциальных уравнений цепи. Дифференциальные уравнения цепи относительно переменных состояния записываются в канонической форме, т.е. представляются решенными относительно первых производных переменных состояния по времени. Число переменных состояния равно порядку системы дифференциальных уравнений исследуемой электрической цепи. -uc - r3·ic + uL = 0; -uc - r3(i1 - iL - i4) + uL = 0; uBx = i1·rBx + uL; i1 =; -uc - r3(- iL - i4) + uL = 0; -uL + i4(r5 + rн) = 0. i4 =. Подставляя значения сопротивлений и упрощая выражение, получим: = – iL10+ uc + uBx/10; uвх = i1rвх + uc + r3ic; -uc - ic·r3 + uL = 0; i4(r5 + rн) = uL; i4 = i1 - ic - iL; -uc - ic·r3 + (i1 - ic - iL)( r5 + rн) = 0; i1 = (uвх - uc - r3ic)/r1 ; Подставляя значения сопротивлений и упрощая выражение, получим: Получаем систему уравнений: = – iL10+ uc + U/10 . ИЛИ = * iL+ *uc + *iL + *uc+ function dUdt=ProcessesRLC(t,f) L=0.55; C=0.000002; U=120*sin(1178*t); dUdt=[(-10/(L*(1+0.1+1/22)))*f(1)+(1/(L*(1+0.1+1/22)))*f(2)+U/10/(L*(1+0.1+1/22));(-220/(10*C+10*C*220/100+220*C))*f(1)+((-1-220/100)/(10*C+10*C*220/100+220*C))*f(2)+U*220/100/(10*C+10*C*220/100+220*C)]; end [t,f]=ode45('ProcessesRLC',[0 0.1],[0 0]); subplot(4,1,1) plot(t,f(:,1),'r-'),axis([0 0.08 -0.2 0.3]),grid subplot(4,1,2) plot(t,f(:,2)),axis([0 0.08 -100 100]),grid %f(:,1)это ток iL согласно системе уравнений Коши %f(:,2)это напряжение на конденсаторе % напряжение на нагрузке ниже!!!!!!!вывод формулы напряжения на нагрузке ниже, после графика subplot(4,1,4) plot(t,120*sin(1178*t),'g'),axis([0 0.08 -150 150]),grid hold on plot(t,0.079*120*sin(1178*t)+0.7936*f(:,2)-7.936*f(:,1)),axis([0 0.08 -200 200]),grid %схема в лекции по матлаб% Моделирование работы двигателя постоянного тока с независимым возбуждением 2пн160Lухл4. к началу 1. Основные соотношения. 2. Механическая характеристика M=[0 10 20 30 40 50 60 70]; w=220/2.4-M*0.8./2.4^2 plot(M,w);axis([0 90 0 100] 3. Пуск двигателя на х.х., мгновенный наброс нагрузки. function w=DPT(t,w) x=[0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5]; y=[0 0 40 40 0 0 0 0 0 0 0]; Mc=interp1(x,y,t); w=((MomDPT(w)-Mc)/0.1) end function M=MomDPT(w) M=(220*2.4-w*2.4^2)/0.8 end [t,w]=ode45(@DPT,[0 2.5],0); plot(t,w) ИЛИ w=dsolve('Dw=660-w*7.2','w(0)=0') ezplot(w,[0 2.5]),axis([0 3 0 100]),grid Моделирование работы асинхронного двигателя АИР63В2 . 1. Основные соотношения. 2. Механическая характеристика Электродвигатель Мощность кВт Об/мин. Ток при 380В, А KПД, % Kоэф. мощн. Iп/ Iн Масса, кг Электродвигатели выпускавшиеся ранее АИР 63 В2 0,55 3000 1,3 75 0,81 5 6,1 4А63В2 4АМ63В2 Электро двигатель Мощность Об/мин.* Ток при 380В, А* KПД, %* Kоэф мощн.* Iп/ Iн Мп/ Мн Мmax/ Мн Момент инерции кгм2* Масса кг* АИР63B2 0,55 кВт 2730 1,3 75 0,81 5 2,2 2,2 0,0009 6,1 λ перегрузочная способность % аир 63в2 sk=0.37 s=0:0.01:1; M=2.*4.23./(sk./s+s./sk)% формула Клосса plot(s,M) gtext('M') gtext('s') или n=0:3000; m=2.* 4.23./(0.37./((3000-n)./3000)+((3000-n)./3000)./0.37) plot(m,n) 3. Пуск двигателя на х.х., мгновенный наброс нагрузки function s = Skoljen(w) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here s=(314-w)./314; end function Mso=Mso(w) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here Mso=2.*0.4725./(0.37./Skoljen(w)+Skoljen(w)./0.37); end function dw=Asinhr(t,w) %UNTITLED4 Summary of this function goes here % Detailed explanation goes here x=[0 2 4 6 8 10]; y=[0 0 0.3 0.3 0 0]; Mc=interp1(x,y,t); dw=(Mso(w)-Mc)./0.0009 end [t,w]=ode45(@Asinhr,[0 10],0);axis([0 10 0 500]); plot(t,w) ИЛИ function dw = dv(t,w) dw=1050/(0.37/(1-w/314)+(1-w/314)/0.37) end [t,w]=ode45(@dv,[0 10],0); plot(t,w) С нагрузкой function dw = dv(t,w) x=[0 2 4 6 8 10]; y=[0 0 0.2 0.2 0 0]; Mc=interp1(x,y,t); dw=1050/(0.37/(1-w/314)+(1-w/314)/0.37)-Mc/0.0009 end [t,w]=ode45(@dv,[0 10],0); plot(t,w) Аппроксимация в Matlab Пример мой: clc u=[0 1 2 3 4 5]; - задан график – зависимость тока от напряжения i=[0 1 16 27 64 125]; subplot(2,1,1); plot(u,i); axis([0 6 0 200]); - строится этот график по заданным таблично точкам p = polyfit(u,i,3) – команда для определения коэффициентов полинома t=0:0.01:10; f=1.2963.*t.^3-2.7937.*t.^2+6.6878.*t-0.8889; - проверка (построение графика по найденным коэффициентам полинома. Графики должны совпасть) subplot(2,1,2) plot(t,f,'r') axis([0 6 0 200]); что выдалось в командном окне: это и есть коэффициенты полинома. Разберём задачу, в которой разрешается использование встроенных  функций. Осуществить аппроксимацию в Matlab табличных данных x = [0, 0.1 , 0.2, 0.3, 0.5] и y = [3, 4.5, 1.7, 0.7, -1] . Применяя метод наименьших квадратов, приблизить ее многочленами 1-ой и 2-ой степени. Для каждого определить величину среднеквадратической ошибки. Построить (на одном листе) графики и заданной таблично функции (ломанная линия) и приближающих ее многочленов 1-ой и 2-ой степени. x = [0, 0.1 , 0.2, 0.3, 0.5]; y = [3, 4.5, 1.7, 0.7, -1]; grid on plot(x, y, '*r'); xi = min(x):0.1:max(x); N = 1; % степень coeff1 = polyfit(x, y, N); y2 = 0; for k=0:N y2 = y2 + coeff1(N-k+1) * xi.^k; end hold on; plot(xi, y2, 'r'); N = 2; coeff2 = polyfit(x, y, N); y3 = 0; for k=0:N y3 = y3 + coeff2(N-k+1) * xi.^k; end hold on; plot(xi, y3, 'g'); std(y-(coeff1(1)*x+coeff1(2))) std(y-(coeff2(1)*x.^2+coeff2(2)*x+coeff2(3))) Вывод: ans = 0.9253 ans = 0.8973 Однако, встречаются задачи, где требуется реализовать аппроксимацию в Matlab без использования специальных функций. Найти у(0.25) путём построения аппроксимирующего полинома методом наименьших квадратов согласно данным: x: 0, 0.1, 0.2, 0.3, 0.5 y: 3, 4.5, 1.7, 0.7, -1 p: 0.5, 0.8, 1.6, 0.8, 0.1 Построить этот полином без учёта весовых коэффициентов с использованием определителя Вандермонда и стандартных операторов. %Задаем массивы данных: x = [0; 0.1; 0.2; 0.3; 0.5]; y = [3;4.5;1.7;0.7;-1]; %Строим матрицу W – матрицу Вандермонда с вырезанным первым столбцом: W = vander(x); W = W(1:5,2:5); %Вычисляем элементы матрицы А как произведение транспонированной матрицы W и просто матрицы W A = W'*W; %Вычисляем элементы вектора b b = W'*y; %Решая систему уравнений Aa = b, находим значения вектор-столбца a: a = inv(A)*b %Это будут коэффициенты аппроксимирующего полинома. %Проверяем, используя методы MATLAB = функцию polyfit: qq = polyfit(x,y,3) %Получаем аппроксимированные значения y: x1 = [-0.2:0.001:0.7]; y1 = a(1)*x1.^3 + a(2)*x1.^2 + a(3)*x1 + a(4); %Строим график функции plot(x,y,'*'); hold on; grid on; plot(x1,y1,'Color','r'); %Находим значение в точке x = 0.25 x2 = 0.25; y2 = a(1)*x2^3 + a(2)*x2^2 + a(3)*x2 + a(4) Вывод: a = 228.1447 -176.0984 22.7745 3.1590 qq = 228.1447 -176.0984 22.7745 3.1590 y2 = 1.4113 Как видите встроенные функции для аппроксимации в Matlab укорачивают алгоритм почти вдвое. Существует также возможность реализации всего алгоритма через одну функцию, но для преподавателей студентов она скорее всего будет не приемлема. С помощью функции lsqcurvefit(fun,x0,xdata,ydata), где: xdata,ydata– табличные значения аппроксимируемой функции; x0 –стартовое значение параметров функции; fun – функция аппроксимации, задаваемая пользователем С аналитически-теоретической стороны, существуют такие виды аппроксимации: • Аппроксимация ортогональными классическими полиномами. • Аппроксимация каноническим полиномом Но на практике их реализацию требуют редко. Вот и вся основная информация по аппроксимации в Matlab, если остались вопросы, задавайте их в комментариях.