Анализ систем автоматизированного управления численными методами

Тема: Аналіз САУ чисельними методами

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

1. Техническое задание на проектирование Задана структурная схема следящей системы постоянного тока

Общий коэффициент передачи разомкнутой следящей системы определяется как:

kv=kиу*kфd*kун*kэму*kдв*kред

Таблица 1.1 Исходные данные для проектирования:

Необходимо: 1. Найти критический коэффициент (kvкрит) передачи следящей системы используя критерий устойчивости Гурвица и метод исключения Гаусса 2. Найти и графически построить переходные функции замкнутой системы для 2х вариантов: а). kv=½kvкрит б). kv=¼kvкрит

используя метод Рунге-Кутта. 3. Вычислить квадратичную интегральную оценку для 2х вариантов по импульсной переходной функции используя метод Симпсона

2. Расчёт критического коэффициента передачи замкнутой следящей системы

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

p (Tкзp+1)(Тупрp+1)[(Tэмp+1)Tяp+1]

тогда характеристический многочлен замкнутой системы равен:

A (p)=p (Ткзp+1)(Tупрp+1)[(Tэмp+1)Tяp+1]+kv=a5p+a4p+a3p+a2p+a1p+a0

a5=Tкз*Тупр*Тэм*Тя;

a4=Тупр*Tэм*Tя+Tкз*Tэм*Tя+Tкз*Тупр*Tя;

a3=Tэм*Tя+Тупр*Tя+Tкз*Tя+Tкз*Тупр;

a2= Tя+Тупр+Tкз; a1=1; a0= kv;

2.1 Методика расчёта критического коэффициента

Критическим коэффициентом kvкрит называется коэффициент передачи разомкнутой системы находящейся на границе устойчивости В соответствии с критерием устойчивости Гурвица составляем определитель Гурвица. Для заданной системы 5-го порядка этот определитель имеет следующий вид:

Система будет находиться на границе устойчивости, если один из диагональных миноров равен 0, а все остальные — положительны. Процедура поиска критического коэффициента передачи состоит в следующем: задавая, последовательно увеличивая значение kv, начиная от kvнач=0.05 с дискретностью kv=0.05 определяем такое его значение при котором хотя бы один из миноров становится равным 0. Следовательно для каждого значения kv нужно вычислить все пять определителей. Вычисление определителей будем проводить методом исключения Гаусса с выбором главного элемента.

2.2 Метод Гаусса с выбором главного элемента

Метод исключения Гаусса состоит из двух основных этапов:

1) прямой ход — основан на приведении матрицы системы к треугольному виду. Это осуществляется последовательным исключением неизвестных из уравнений системы. 2) обратный ход — с помощью найденной треугольной матрицы последователь, но вычисляем искомые неизвестные. Рассмотрим применение метода Гаусса для системы 3-го порядка (прямой ход):

а11×1+а12×2+а13×3=b1, (1)

а21×1+а22×2+а23×3=b2, (2)

а31×1+а32×2+а33×3=b3. (3)

Домножим уравнение (1) на (- а21/а11) и прибавим его к уравнению (2). Затем, умножив уравнение (1) на (- а31/а11) и прибавив результат к уравнению (3), получим равносильную систему уравнений вида а11×1+а12×2+а13×3=b1, (4) а22×2+а23×3=b2, (5)

а32×2+а33×3=b3, (6)

где

а22 = а22 — (а21/а11)* а12;

b2 = b2 — (а21/а11)* b1;

а23 = а23 — (а21/а11)* а13;

b3 = b3 — (а31/а11)* b1;

а32 = а32 — (а31/а11)* а12;

а33 = а33 — (а31/а11)* а13;

Домножим уравнение (5) на (- а32/а22) и складываем его с уравнением (6)

а11×1+а12×2+а13×3=b1, (7) а22×2+а23×3=b2, (8) а33×3=b3, (9)

где а33 = а33 — (а32/а22)* а23;

b3 = b3 — (а32/а22)* b2;

Мы привели матрицу системы к треугольному виду. На этом заканчивается прямой ход метода Гаусса. Обратный ход начинается с решения (9) уравнения системы:

х3 = b3/a33;

Используя это значение, можно найти х2 из (8) уравнения, а затем х1 из (7):

х2 = 1/a22*(b2 — a23×3), х1 = 1/a11*(b1 — a12×2 — a13×3);

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

3. Текст программы

program Krit;

uses crt;

const n=5; Tem=0.131; Tya=0.05; Typr=0.076; Tkz=0.148;

type massiv= array [1.n, 1. n] of real;

var a, b: massiv; Kv, diag: real; i, j: integer;

procedure Vivod;

begin textcolor (9); writeln (' Коэффициенты определителя Гурвица :'); textcolor (7); writeln (' a1=', 1); writeln (' a2=', Tya+Typr+Tkz:9:8); writeln (' a3=', Tem*Tya+Typr*Tya+Tkz*Tya+Tkz*Typr:9:8); writeln (' a4=', Typr*Tem*Tya+Tkz*Tem*Tya+Tkz*Typr*Tya:9:8); writeln (' a5=', Tkz*Typr*Tem*Tya:9:8);

end;

Обнуление элементов стоящих под главной диагональю

procedure Obnylenie (m:integer; a: massiv; var b: massiv);

var i, j: integer;

begin

clrscr; for i:=m+1 to n do for j:=m+1 to n do b[i, j]: =a[i, j]-a[i, m]/a[m, m]*a[m, j];

end;

Вычисление определителей

procedure Diag_opred (b:massiv; var sum: real);

var i, j: integer;

begin clrscr; Vivod; writeln; textcolor (9); writeln (' Вычисленные определители:'); sum:=1; for i:=1 to n do for j:=1 to n do if i=j then begin sum:=sum*b[i, i]; textcolor (7); writeln (' ', sum:7:9); end;

end;

Формирование определителя Гурвица

procedure Matrix (Kv:real;var a: massiv);

begin clrscr; for i:=1 to n-3 do begin a[2*i, i+2]: =1; a[2*i, i+1]: =Tem*Tya+Typr*Tya+Tkz*Tya+Tkz*Typr; a[2*i, i]: =Tkz*Typr*Tem*Tya; end; for i:=1 to n-2 do begin a[2*i-1,i+2]: =Kv; a[2*i-1,i+1]: =Tya+Typr+Tkz; a[2*i-1,i]: =Typr*Tem*Tya+Tkz*Tem*Tya+Tkz*Typr*Tya; end;

end;

Основное тело программы

begin Kv:=0; repeat begin Kv:=Kv+0.05; Matrix (Kv, a); end; b:=a; for i:=1 to n-1 do Obnylenie (i, b, b); begin Diag_opred (b, diag); end; until diag<=0; begin writeln; textcolor (9); writeln (' Критический коэффициент равен:'); textcolor (7); writeln (' Kv =', Kv:5:2); readkey; end;

end.

3.1 Описание переменных

В программе использованы такие переменные: a-массив чисел, из элементов которого формируем определитель Гурвица; n-размерность массива; i-номер строки; j-номер столбца; b-массив чисел. При помощи его, изменяя номера строк и столбцов приводим матрицу к треугольному виду; Tem-электро-механическая постоянная времени двигателя; Tyaпостоянная якорной цепи двигателя; Typr-постоянная времени обмотки управления ЭМУ; Tkz-постоянная времени короткозамкнутой обмотки ЭМУ; m-некоторая переменная принимающая значение от 1 до 5; sum-счетчик, используемый для вычисления диагональных определителей; Kv-дискретность или шаг; diag-один из вычисленных диагональных определителей;

3.2 Описание программы

В процедуре «Vivod» осуществляется расчет коэффициентов a1, a2,a3,a4,a5 определителя Гурвица и вывод значений на экран; В процедуре «Obnylenie» задавая некоторой переменной «m» значения от 1 до 5 мы в цикле изменяя номера строк и столбцов приводим матрицу размерностью 55 к треугольному виду, используя дополнительный массив b. В этом массиве накапливается значение элемента i-ой строки и j-го столбца; В процедуре «Diag_opred» происходит вычисление диагональных определителей. Присваивая некоторому счетчику «sum» значение равное 1, мы при первом выполнении цикла получаем 1-ый главный диагональный минор определителя Гурвица, следующий получаем умножая 1-ый на 2-ой и так далее; В процедуре «Matrix» — выстраиваем сам определитель Гурвица. Для облегчения в цикле зададим значение только переменной i; Tеперь рассмотрим основное тело программы: пусть первоначально Kv=0, затем выполняем вычисление диагональных определителей, каждый раз увеличивая шаг на 0.05 т. е. Kv=Kv+0.05 до тех пор пока один из них не станет меньше ноля (diag<0). Как только это условие выполнится выводим на печать значение cоответствующего Kv-это и будет критический коэффициент передачи замкнутой следящей системы.

3.3 Результаты работы программы

переходный функция импульсный квадратичный Коэффициенты определителя Гурвица: a1=1;

a2=0.27 500 000;

a3=0.1 710 700;

a4=0.76 100;

a5=0.2 110;

Вычисленные определители: 0.760 998; 0.7 217; 0.1 637; -0.2; -0.32;

Критический коэффициент равен: kv=14,40;

4. Определение переходной функции h (t) следящей системы и её показатели качества

4.1 Методы нахождения передаточной функции

Передаточная функция следящей системы равна:

kv

k (p)=

p (Tкзp+1)(Тупрp+1)[(Tэмp+1)Tяp+1]

где

A (p)=p (Ткзp+1)(Tупрp+1)[(Tэмp+1)Tяp+1]+kv=a5p+a4p+a3p+a2p+a1p+a0

характеристический многочлен разомкнутой системы Т.к. а0=b0=kv мы можем переписать данную функцию в таком виде:

b0 Y (p)

k (p)= a5p+a4p+a3p+a2p+a1p+a0 X (p)

На основании данной функции запишем дифференциальное уравнение замкнутой системы. Оно имеет вид:

a5*dy/dt +a4*dy/dt +a3*dy/dt +a2*dy/dt +a1*d y/dt +a0*y=b0*x;

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

dy/dt +a4*dy/dt +a3*dy/dt +a2*dy/dt +a1*dy/dt +a0*y=b0*x;

где ai =ai/a5,(i=0,1,…, 4); a0=b0=b0/a5;

Введём новую переменную z, тогда z1=y, z2=y, z3=y, z4=y, z5=y.

Преобразовав выражение представим систему в таком виде: dz1/dt=z2, dz2/dt=z3, dz3/dt=z4, dz4/dt=z5, dz5/dt=b0x-a0z1-a1z2-a2z3-a3z4-a4z5. В нашем случае используется метод Э2, начальные условия:

z1(0)=z2(0)=z3(0)=z4(0)=z5(0). В результате решения необходимо получить таблицу значений для функции z1(t)(при двух значениях: kvкрит/2 и kvкрит/4) h (t)=z1(t) — переходная функция — реакция системы на единичное ступенчатое воздействие, причем её значения должны затухать вокруг 1. Время t, при котором данная функция попадает в 5% коридор отклонения от значения 1 и не выходит из него называется tрег-время регулирования. По данным таблиц строим графики z1(t)=h (t). По найденным далее значениям t и t, соответствующих kv крит/2 и kv крит/4 вычисляем интегральные квадратичные оценки:

I=

I=

Определение значения интегралов производим по методу Симпсона.

4.2 Метод Рунге-Кутта

Наиболее распространенным среди одношаговых методов является метод Рунге-Кутта. На его основе могут быть построены разностные схемы разного порядка точности. Приведем схему метода четвертого порядка.

Таким образом, метод Рунге-Кутта требует на каждом шаге четырехкратного вычисления правой части уравнения f (x, y), что приводит к большему объёму вычислений, однако это окупается повышенной точностью, что дает возможность проводить счет с большим шагом.

5. Вычисление интергальной квадратичной оценки по импульсной переходной функции методом Симпсона

Для наглядности рассмотрим в чем заключается метод Симпсона. Для этого разобьем отрезок интегрирования [a, b] на четное число n равных частей с шагом h. На каждом отрезке [x0,x2],[x2,x4],…,[xi-1,xi+1],…,[xn-2,xn] подинтегральную функцию (х) заменим интерполяционным многочленом второй степени:

(x) i (x)=aix+bix+ci, xi+1 x xi-1

В качестве i (x) можно принять интерполяционный многочлен Лагранжа второй степени, проходящий через точки Mi-1(xi-1,yi-1), Mi (xi, yi),

Mi+1(xi+1,yi+1): (x-xi)(x-xi+1) (x-xi-1)(x-xi+1) (x-xi-1)(x-xi)

i (x)= yi-1+ yi+ yi+1; (xi-1-xi)(xi-1-xi+1) (xi-xi-1)(xi-xi+1) (xi+1-xi-1)(xi+1-xi)

Проведя вычисления для каждого элементарного отрезка [xi-1,xi+1], просуммируем полученные выражения:

S=h/3(y0+4y1+2y2+4y3+2y4+…+2yn-2+4yn-1+yn).

Данное выражение для S принимается в качестве значения определенного интеграла:

(x)dxh/3[y0+4(y1+y3+…+yn-1)+2(y2+y4+…+yn-2)+yn].

Полученное соотношение и есть формула Симпсона. Её также можно получить и другими способами, например, комбинированием формул прямоугольников и трапеций или двукратным применением метода трапецій при разбиении отрезка [a, b] на части с шагами h и 2h. Главный член погрешности метода Симпсона имеет вид:

Rn=(-h/180) (x).

6. Текст программы

program kurs2;

uses crt, dos;

const Kv=7.2; b0=7.2; a0=7.2; a1=1; a2=0.275; a3=0.17 107; a4=0.761; a5=0.211; h=0.05;

type mas=array [1.5,0.100] of real;

var

integral, d, c0,c1,c2,c3,c4:real;

k, n, j:integer;

Z, K1, K2,K3,K4:mas;

Y:array [1.3] of real;

c:char;

begin

clrscr;

integral:=0;

d:=b0/a5;

c0:=a0/a5;

c1:=a1/a5;

c2:=a2/a5;

c3:=a3/a5;

c4:=a4/a5;

for j:=1 to 5 do

Z[j, 0]: =0;

for j:=1 to 5 do begin K1[j, 0]: =0; K2[j, 0]: =0; K3[j, 0]: =0; K4[j, 0]: =0; end; begin n:=0; repeat for j:=1 to 4 do K1[j, n]: =h*Z[j+1,n]; K1[5,n]: =h*(d-c0*Z[1,n]-c1*Z[2,n]-c2*Z[3,n]-c3*Z[4,n]-c4*Z[5,n]);

for j:=1 to 4 do K2[j, n]: =h*(Z[j+1,n]+K1[j+1,n]/2); K2[5,n]: =h*(d-c0*z[1,n]-c1*z[2,n]-c2*Z[3,n]-c3*Z[4,n]-c4*Z[5,n]); for j:=1 to 4 do K3[j, n]: =h*(Z[j+1,n]+K2[j+1,n]/2); K3[5,n]: =h*(d-c0*Z[1,n]-c1*Z[2,n]-c2*Z[3,n]-c3*Z[4,n]-c4*Z[5,n]); for j:=1 to 4 do K4[j, n]: =h*(Z[j+1,n]+K3[j+1,n]); K4[5,n]: =h*(d-c0*Z[1,n]-c1*Z[2,n]-c2*Z[3,n]-c3*Z[4,n]-c4*Z[5,n]); for j:=1 to 5 do begin Z[j, n+1]: =Z[j, n]+(K1[j, n]+2*K2[j, n]+2*K3[j, n]+K4[j, n])/6; end; y[k]: =sqr (z[2,n]); if k=3 then begin integral:=integral+(h/3)*(y[1]+4*y[2]+y[3]); y[1]: =y[3]; k:=1; end; textcolor (3); writeln ('No итерациии — ', n,' Z1 = ', Z[1,n]: 3:8,' T= ',((n*h)-0.05):3:2,' c.', ' Интеграл = ', integral:3:8); normvideo; c:=readkey; k:=k+1; n:=n+1; until c='e' end;

end.

6.1 Описание переменных

В программе использованы такие переменные:

Массив Z — решения диф. уравнений.

Массив К1, К2, К3, К4 — коэффиенты для решения диф. уравнений по методу Рунге-Кутта.

С0,с1,с2, с3, d — откорректированные коэффициенты для решения диф. уравнения 5-го порядка.

У — массив для нахождения интегральной оценки.

J, n, k — счетчики столбцов и строк массивов.

Integral — переменная, в которую на каждом шаге заносится значение интегральной оценки.

6.2 Описание работы программы

На первом шаге программа обнуляет все массивы, вычисляет коэффициенты. Далее программа в цикле производит решение уравнений по методу Рунге-Кутта, при этом коэффициенты К на каждом шаге находятся сразу для всех уравнений.

Далее программа находит интегральную оценку и выводит результаты на экран. Итеррации будут продолжаться до тех пор пока мы не нажмём клавишу `e'.

Шаг на каждой итеррации увеличивается на 0,05. Интеграл считается на каждом 3-м шаке.

6.3 Результаты работы программы

Таблица 8.1

Выводы

1. Для следящей схемы показанной на рисунке 1.1, был расчитан критический коэффициент передачи, который равен 14.40.

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

3. Для нахождения переходной характеристики необходимо решить систему дифференциальных уравнений, используя метод Рунге-Кутта, имеющий четвёртый порядок точности. Другими словами, метод Рунге-Кутта даёт при незначительных затратах на вычисления более точные данные чем другие методы (метод Эйлера). Решение задачи показало целесообразность применения данного метода.

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

5. При построении переходных характеристик для значений kv/2 и kv/4, мы определили, что при уменьшении kv система становится устойчивой. При значении kv/4 колебания затухают быстрее, нежели при kv/2

Список используемых источников

1. Турчак Л. И. Основы численных методов. — М.: Высшая школа, 1987 г.

2. Заворыкин В. М. Численные методы. — М.: Высшая школа, 1991 г.

3. Виньяминов Б. Б., Рогачёв А. И. Методы приближённых исследований систем САУ на ЭВМ, 1989 г.