Козлов Олег Степанович
Кондаков Дмитрий Евгеньевич
Скворцов Леонид Маркович
Тимофеев Константин Александрович
Ходаковский Виктор Владимирович

Программный комплекс
«Моделирование в технических устройствах» (ПК «МВТУ»)

Аннотация: Рассматриваются основные возможности программного комплекса «МВТУ», предназначенного для исследования динамики и проектирования технических систем и устройств различного назначения. Математические модели исследуемых систем формируются в виде структурных схем, элементы которых описываются входо-выходными соотношениями, дифференциально-алгебраическими и разностными уравнениями, логическими условиями. Для исследования и проектирования таких систем в ПК «МВТУ» реализованы методы имитационного моделирования, поисковой оптимизации и статистического анализа, а также частотные и корневые методы анализа и синтеза динамических систем.

  1. Введение
  2. Построение математических моделей
    1. Формирование моделей из типовых блоков
    2. Язык программирования
    3. Реализация дифференциально-алгебраических уравнений
  3. Моделирование
  4. Оптимизация
  5. Анализ
  6. Синтез
    1. Частотный метод
    2. Корневой метод
  7. Исследование адаптивного ПИ-регулятора
  8. Заключение
  9. Литература и Интернет

Введение

Программный комплекс «МВТУ» [1–5] предназначен для исследования динамики и проектирования самых разнообразных систем и устройств. По своим возможностям он является альтернативой аналогичным зарубежным программным продуктам Simulink, VisSim и др. Удобный редактор структурных схем, обширная библиотека типовых блоков и встроенный язык программирования позволяют реализовывать модели практически любой степени сложности, обеспечивая при этом наглядность их представления. ПК «МВТУ» успешно применяется для проектирования систем автоматического управления, следящих приводов и роботов-манипуляторов, ядерных и тепловых энергетических установок, а также для решения нестационарных краевых задач (теплопроводность, гидродинамика и др.). Широко используется в учебном процессе, позволяя моделировать различные явления в физике, электротехнике, в динамике машин и механизмов, в астрономии и т.д. Может функционировать в многокомпьютерных моделирующих комплексах, в том числе и в режиме удаленного доступа к технологическим и информационным ресурсам.

ПК «МВТУ» реализует следующие режимы работы:

К достоинствам ПК «МВТУ» относятся:

Для отечественных пользователей удобство работы с ПК «МВТУ» обусловлено также русскоязычным интерфейсом и наличием обширной документации на русском языке. Учебная и демонстрационная версии ПК «МВТУ» вместе с полной документацией и набором демонстрационных примеров распространяются свободно [1]. В учебной версии есть ограничения на сложность модели: порядок дифференциальных уравнений не выше 30, а число блоков не более 100. В демонстрационной версии таких ограничений нет, но модель нельзя сохранить.

Рассмотрим основные возможности ПК «МВТУ» на конкретных примерах, которые есть среди демонстрационных примеров в каталоге Demo.

1. Построение математических моделей

1.1. Формирование моделей из типовых блоков

На рис. 1 показаны основные окна ПК «МВТУ» на примере моделирования динамики САР ядерного реактора с релейным регулятором (файл \Demo\Лабораторные_работы\reac_lab.mrj). В верхней части рис. 1 представлено главное окно, содержащее командное меню, панель инструментов и «линейку» типовых блоков, а в нижней части – схемное и графическое окна. Сформированная в схемном окне математическая модель динамики САР описывает нейтронно-физические процессы в ядерном реакторе с использованием точечной модели кинетики (реализованной типовым блоком Кинетика нейтронов) с учетом шести групп ядер-предшественников запаздывающих нейтронов.

Рис. 1. Основные окна ПК «МВТУ».

Для построения математических моделей в ПК «МВТУ» используются типовые элементы – блоки. Каждый типовой блок реализует математическую модель того или иного явления, процесса или устройства. Связи блока с другими блоками структурной схемы задаются через его входы и выходы. Параметры блока определяют коэффициенты в уравнениях математической модели. Библиотека типовых блоков состоит из Общетехнической и ряда специализированных библиотек, доступ к которым осуществляется из «линейки» типовых блоков выбором закладки с соответствующим названием.

Общетехническая библиотека содержит 165 типовых блоков, сгруппированных в 11 каталогов (Источники, Данные, Операции математические, Векторные операции, Субструктуры, Динамические звенья, Нелинейные звенья, Логические звенья, Функции математические, Ключи, Дискретные звенья).

В ПК «МВТУ» включены также 13 специализированных библиотек (Кинетика нейтронов, Свойства (вода и пар, аммиак, смеси идеальных газов, углеводородные газы с учетом сжимаемости), Статистика, Внешние, Анимация, Контроль и управление, Теплопроводность, Реакторные блоки, Гидроавтоматика, Электромашины, Логика АСУ ТП, Элементы ПХГ (подземное хранилище газа), Роботы).

Ряд моделей, сформированных в процессе отладки и эксплуатации ПК «МВТУ», сохранен в отдельных каталогах в виде субмоделей (макроблоков). Среди них – модели приводных устройств и элементов электрических схем. Фактически эти каталоги – дополнительные специализированные библиотеки, которые также могут быть использованы для формирования математических моделей динамических систем.

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

Формирование математической модели предполагает задание параметров каждого блока структурной схемы. Часто параметры функционально зависят друг от друга и от некоторых общих для всей модели параметров. Для удобства задания параметров в ПК «МВТУ» они подразделяются по области действия на локальные и глобальные.

Локальные параметры являются численными характеристиками элементарного блока. Сфера их действия ограничена математической моделью этого блока.

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

В ПК «МВТУ» с помощью типового блока В память могут быть заданы глобальные переменные модели. В качестве параметра блока В память указывается имя, которое дается входной переменной этого блока. Глобальные переменные задаются, а затем используются в любой части структурной схемы. Типовой блок Из памяти позволяет выбрать одну из глобальных переменных, которая становится выходом этого блока.

Глобальные переменные широко используются в различных режимах работы программного комплекса. Во-первых, такой механизм является удобным средством для проведения невидимых линий связи между любыми уровнями структурной схемы, что позволяет улучшить читаемость сложных многосвязных моделей. Во-вторых, глобальные переменные используются для задания критериев качества в режиме работы ОПТИМИЗАЦИЯ. В‑третьих, эти переменные служат для задания входов и выходов в режимах АНАЛИЗ и СИНТЕЗ. Наконец, они используются при обмене данными с виртуальными мнемосхемами и пультами управления в режиме КОНТРОЛЬ И УПРАВЛЕНИЕ.

1.2. Язык программирования

Учитывая, что невозможно сформировать абсолютно полную библиотеку типовых блоков, в ПК «МВТУ» включены средства, позволяющие пользователю расширить состав личной библиотеки путем создания новых типов блоков. К таким средствам относятся типовые блоки Новый и Язык программирования. Блок Язык программирования позволяет создавать модели с помощью алгоритмического языка, аналогичного языку системы MATLAB. Диалоговое окно этого блока – окно текстового редактора алгоритмов, в котором пользователь записывает математическую модель в виде, близком к естественной записи. Язык позволяет задавать алгебраические соотношения и дифференциальные уравнения, а также выполнять операции с действительными и комплексными матрицами и векторами, логическими переменными, геометрическими точками, полиномами. Возможны и более сложные типы данных – записи и ссылки. Среди конструкций языка – условные и безусловные переходы, циклы, пользовательские функции и процедуры. Интерпретатор языка предварительно компилирует программу во внутренний псевдокод, что ускоряет выполнение задачи.

В качестве примера рассмотрим построение модели одномерной теплопроводности, описываемой уравнением

при заданных граничных и начальных условиях (файл \Demo\Язык_программирования\ Теплопроводность.mrj). Применяя для дискретизации этого уравнения метод прямых, получим систему ОДУ

Текст алгоритма, реализующего решение рассматриваемой задачи, записан в окне редактора блока Язык программирования в следующем виде:

Листинг 1

input u,v;   //входы блока (граничные условия)
init T=N#20; //начальные условия
const alfa=0.02, K=alfa*(N-1)^2,
x=linspace(0,1,N); 
//x - массив узлов равномерной сетки
T[1]=u; T[N]=v;
for (i=2,N-1) T'[i]=K*(T[i-1]-2*T[i]+T[i+1]);
output x[N],T[N]; //векторные выходы блока

(комментарии следуют за двойным слешем — "//"). Число узлов N задается в окне глобальных параметров. Для отображения распределения температуры используется типовой блок График Y от X, на входы которого подаются векторные сигналы x и T.

Рассмотрим еще один пример – фильтрацию сигналов с использованием быстрого преобразования Фурье (БПФ) (файл \Demo\Язык_программирования\Фильтр_БПФ.mrj). Пусть исходный сигнал – сумма двух гармоник, а на него накладывается нормальный белый шум, сопоставимый по величине с исходным сигналом. Тогда алгоритм фильтра БПФ, записанный с помощью языка программирования, можно представить в виде

Листинг 2

t=linspace(0,1023,1024)/4;
x=sin(5*pi*t/128)+0.9*cos(pi*t/8);
xn=1024#0;
for (k=1,1024) xn[k]=x[k]+2*randg(0,1);
// t - массив значений времени
// x - исходный сигнал
// xn - зашумленный сигнал
y=fft(xn);  //прямое БПФ
py=abs(y).*abs(y);
Porog=4e4;
// py - спектральная плотность сигнала
// Porog - порог для спектральной плотности
pz=1024#(0,0);
for (k=1,1024)
if py[k]<Porog then pz[k]=0+0i else pz[k]=y[k];
z=real(ifft(pz)); //обратное БПФ
// z - фильтрованный сигнал
f=linspace(1,511,511)/256;
// массив значений частоты
for (k=1,511) E[k]=py[k+1]*1e-4;
E_Porog=511#(Porog*1e-4);
// нормированные значения спектр. плотности и порога
output t[1024],xn[1024],x[1024],f[511],
       E[511],E_Porog[511],z[1024];

На выходе блока получаем массивы значений сигналов, а также массив значений спектральной плотности. Результаты, отображенные с помощью блоков График Y от X, показаны на рис. 2.

Рис. 2. Фильтрация сигнала с помощью БПФ: а - зашумленный сигнал (зеленая линия) и исходный сигнал (красная линия); б –спектральная плотность (красная линия) и порог для спектральной плотности (синяя линяя); в – исходный сигнал (красная линия) и фильтрованный сигнал(синяя линия).

Как правило, сложная техническая система является гибридной системой, поведение которой описывается не только дифференциально-алгебраическими уравнениями, но и логическими условиями, задающими переходы из одного дискретного состояния в другое. Для моделирования таких систем применяют специальные программные средства, такие как Simulink+Stateflow или MVS [6]. Язык программирования и большой набор дискретных, переключательных и логических элементов позволяют реализовывать гибридные модели также и в ПК «МВТУ».

Рассмотрим пример, включенный в состав демонстрационных примеров пакета Stateflow под названием Stick-Slip Friction Demonstration. Моделируется движение бруска, прикрепленного к пружине, под действием внешней силы и с учетом сухого трения (такая модель подробно рассмотрена в [7]). Эта же модель, реализованная в ПК «МВТУ» (файл \Demo\Язык_программирования\ Сухое_трение_анимация.mrj), показана на рис. 3. В этом примере моделирование выполняется с использованием типового блока Анимация, демонстрирующего движение бруска в реальном масштабе времени.

Рис. 3. Модель движения бруска с сухим трением.

Параметры модели: M – масса бруска; K – упругость пружины; Fstatic – сила, которую нужно приложить, чтобы сдвинуть брусок с места; Fsliding – сила трения скольжения (Fstatic ³Fsliding). Изменение состояния бруска описывается следующими логическими условиями: если брусок остановился или был в состоянии покоя и при этом суммарная сила Fsum не превышает по модулю Fstatic, то брусок остается в состоянии покоя; в противном случае он находится в состоянии движения. В соответствии с этими условиями алгоритм вычисления скорости бруска, реализованный в блоке Язык программирования, имеет вид

Листинг 3

input Fsum;
init V=0;
var Vold=V;
if (V*Vold<=0) and (abs(Fsum)<=Fstatic) 
then begin V'=0; V=0 end           //покой
else V'=(Fsum-Fsliding*sign(V))/M; //движение
if goodstep then Vold=V;     
//Vold - скорость на предыдущем успешном шаге
output V;

Здесь goodstep – системная логическая переменная, которая принимает значение true в случае успешного шага интегрирования и false в противном случае. На каждом шаге интегрирования производится оценка ошибки. Шаг считается успешным (goodstep=true), если эта оценка меньше допустимой. В противном случае (goodstep=false) происходит возврат к предыдущему шагу, после чего выполняется шаг меньшего размера. Изменение состояния (покой или движение) фиксируется только при успешном шаге. Благодаря использованию языка программирования модель оказалась более простой и наглядной, чем аналогичная модель в среде Simulink+Stateflow. Отметим, что и время счета такой модели в ПК «МВТУ» в несколько десятков раз меньше, чем в Simulink (для определения времени счета следует использовать пример без анимации – Сухое_трение.mrj).

Язык программирования ПК «МВТУ» позволяет каждому пользователю создавать собственные функции на основе стандартных конструкций языка. Такие функции могут быть объединены в библиотеки, вызов которых осуществляется оператором вида

Листинг 4

include "filename.txt";

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

Листинг 5

include "procedures.txt";
input B[M],A[N]; 
//массивы коэффициентов передаточной функции
w=logspace(-1,3,200); 
//массив частот в логарифмическом масштабе
Nyquist(B,A,w,X,Y);
output X[200],Y[200]; //массивы точек годографа

где процедура расчета точек годографа содержится в файле procedures.txt и имеет вид

Листинг 6

procedure Nyquist(B[1],A[1],w[1],out Re[1],out Im[1])
  var z:complex,i:integer;
  for (i=1,cols(w))
  begin
    z=1j*w[i]; z=polyval(B,z)/polyval(A,z);
    Re[i]=real(z); Im[i]=imag(z);
  end;
end;

Размерности массивов B и A (M и N) задаются в окне глобальных параметров модели. Для отображения годографа следует использовать типовой блок График Y от X, на вход которого подаются векторные сигналы X и Y. Расчет производится на каждом шаге интегрирования, что позволяет наблюдать изменение годографа в процессе моделирования нестационарных и нелинейных систем.

1.3. Реализация дифференциально-алгебраических уравнений

При использовании ПК «МВТУ», а также аналогичных программных комплексов (Simulink, VisSim и др.) наиболее просто реализуются модели, элементы которых описываются соотношениями «вход-выход» либо «вход-состояние-выход», что соответствует явным алгебраическим зависимостям и дифференциальным уравнениям в форме Коши. Однако часто модели технических устройств описываются системами дифференциально-алгебраических уравнений (ДАУ), приведение которых к форме Коши затруднено или невозможно. Поэтому в состав библиотеки ПК «МВТУ» включен блок F(Y)=0 (решение алгебраических уравнений). Блок векторизован (вход и выход могут быть векторами), а его выход рассчитывается исходя из условия равенства нулю входного сигнала.

Система ДАУ может быть представлена в полуявной форме

(1.1)

x' = f(x, y, t), 0 = g(x, y, t)

либо в неявной форме

(1.2)

F(x', x, y, t) = 0

Схемы реализации уравнений вида (1.1) и (1.2) в ПК «МВТУ» показаны на рис. 4, где предполагается, что все линии связи – векторные. Здесь для формирования функций f и g используется блок Язык программирования, а для формирования функции F – блок Новый. При реализации уравнения (1.2) используется также блок Демультиплексор, который «расщепляет» вектор неизвестных размерности n + m на вектор x' размерности n и вектор y размерности m.

Рис. 4. Схемы реализации систем ДАУ: а – в полуявной форме (1.1); б – в неявной форме (1.2).

Таким образом, в ПК «МВТУ» можно реализовать практически любую модель, описываемую системой ДАУ. Например, можно построить модель электрической схемы, придав ей вид, показанный на рис. 5 (файл \Demo\Электротехника\Выпрямитель.mrj). Все элементы этой схемы являются замаскированными макроблоками, соединенными между собой векторными линиями связи. Векторные переменные имеют две компоненты: напряжение и ток. Модели некоторых элементов схемы показаны на рис. 6.

Рис. 5. Пример построения модели электрической схемы: а – схема выпрямителя; б – эквивалентная схема трансформатора (макроблок).

Рис. 6. Модели элементов электрической схемы (синяя линия – напряжение, красная – ток).

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

Приведенный пример показывает, что с помощью макроблоков и блоков В память, Из памяти и F(Y)=0 в ПК «МВТУ» достаточно просто реализуется контурный метод расчета электрических схем. Реализация некоторых других физических моделей затрудняется тем, что линии связи между блоками могут быть только направленными. В последующих версиях ПК «МВТУ» предполагается ввести также и ненаправленные (бинаправленные) линии связи, что, в частности, удобно для реализации узлового метода расчета электрических схем.

2. Моделирование

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

Для решения дифференциальных уравнений в ПК «МВТУ» реализованы 10 явных и 6 неявных методов, среди которых есть новые оригинальные методы [8, 9]. Известно, что классические явные методы (Рунге-Кутты, Адамса и др.) неэффективны при интегрировании жестких систем, поэтому для решения таких задач обычно используют неявные методы. В ПК «МВТУ» реализованы новые явные адаптивные методы [8], параметры которых автоматически настраиваются на решаемую задачу, благодаря чему они позволяют эффективно решать многие жесткие системы. Кроме этого, один из адаптивных методов (Адаптивный 5) обеспечивает точное воспроизведение быстро осциллирующих решений при наличии собственных значений якобиана вблизи мнимой оси. Результаты решения множества тестовых и прикладных задач показали, что реализованные в ПК «МВТУ» методы позволяют быстро и качественно правильно решать задачи разных типов (жесткие, осциллирующие, локально неустойчивые, разрывные). Некоторые результаты тестовых испытаний приведены в [10].

ПК «МВТУ» позволяет также эффективно решать системы ДАУ. При использовании явных методов алгебраическая подсистема решается независимо от дифференциальной, при этом можно применять один из трех методов (простых итераций, Ньютона-Рафсона, Бройдена). При использовании неявных методов алгебраическая и дифференциальная подсистемы решаются совместно, что позволяет решать системы ДАУ высших индексов. Система ДАУ (1.1) имеет индекс 1, если матрица частных производных ∂g(x, y, t) / ∂y обратима в любой точке на траектории решения. В этом случае можно аналитически либо численно исключить из уравнений вектор y, приведя таким образом систему к форме Коши. Если матрица ∂g / ∂y вырождена, то приведение к форме Коши невозможно, такие системы имеют индекс 2 и выше и называются системами ДАУ высших индексов (подробнее об индексе ДАУ см. в [11]). Простейшая схема, описываемая системой ДАУ индекса 2, показана на рис. 7. Эта схема осуществляет дифференцирование входного сигнала, обеспечивая более высокую точность, чем типовой блок Дифференцирование. Подобные схемы позволяют строить обратные модели (данная схема является обратной моделью для интегратора).

Рис. 7. Модель, описываемая системой ДАУ индекса 2 (дифференцирование входного сигнала).

Системы ДАУ высших индексов часто возникают при решении задач механики, теории управления, электротехники и т.д. [11]. Например, электрическая схема, представленная на рис. 5, описывается системой ДАУ индекса 2. Однако Simulink, Vissim и многие другие известные системы моделирования позволяют решать ДАУ только индекса 1. ПК «МВТУ» имеет очевидное преимущество перед этими ПК, позволяя решать системы ДАУ высших индексов.

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

Удобным средством просмотра результатов в процессе моделирования является блок Параметры горячей линии, который инициализируется после 2-х кратного щелчка левой клавишей мыши по линии связи. В результате появляется окно с текущими значениями сигналов (несколькими для векторной линии). Это позволяет определить значение любой переменной в любой момент модельного времени, что упрощает отладку сложных моделей. В ПК «МВТУ» есть и другие средства отладки. Например, можно найти все алгебраические контуры с указанием блоков, удаление которых размыкает эти контуры. Можно включить режим отладки, тогда после каждого шага интегрирования можно посмотреть значения всех переменных. Средства диагностики обеспечивают выдачу диагностических сообщений («Переполнение», «Деление на 0» и т.д.) с указанием конкретного блока, где произошла данная ошибка.

Для оперативного управления процессом моделирования используется режим КОНТРОЛЬ И УПРАВЛЕНИЕ, позволяющий с помощью специальной библиотеки создать Панель управления с расположенными на ней виртуальными аналогами переключателей, ручных регуляторов, лампочек и измерительных приборов. Все эти устройства можно связать с конкретными переменными модели. Пример такой панели управления показан на рис. 8.


 

Рис. 8. Панель управления.

3. Оптимизация

При решении задачи оптимизации в ПК «МВТУ» необходимо задать оптимизируемые параметры (обозначим их p1, ..., pn) и наложить на них ограничения вида

pi min ≤ pi ≤ pi max , i = 1, ..., n.

Эти параметры задаются как глобальные. Критерии качества вычисляются в процессе моделирования. Они формируются типовыми блоками и задаются как глобальные переменные с помощью блоков В память. На критерии (обозначим их Q1, ..., Qm) также накладываются ограничения вида

Qi min ≤ Qi ≤ Qi max , i = 1, ..., m.

Задача оптимизации формулируется следующим образом: найти значения оптимизируемых параметров, при которых выполняются все заданные ограничения. Для решения этой задачи на основе нормированных частных критериев программа оптимизации формирует общий критерий, который минимизируется с помощью поисковых алгоритмов (некоторые из них приведены в [12]). Оптимизации заканчивается, если выполняется хотя бы одно из трех условий: 1) выполняются все ограничения; 2) приращения по параметрам меньше заданных; 3) превышено максимальное число моделирований. Если на какой либо из критериев наложить заведомо невыполнимые ограничения, то в зависимости от этих ограничений он будет минимизирован либо максимизирован. Подобным образом можно решать задачи условной оптимизации, а в общем случае – задачи многокритериальной оптимизации.

В качестве примера рассмотрим оптимизацию параметров ПИД-регулятора с передаточной функцией

для системы с объектом, имеющим передаточную функцию

.

(файл \Demo\Оптимизация\Opt_PID.mrj). Структурная схема, подготовленная для решения этой задачи в ПК "МВТУ" показана на рис. 9, где объект представлен макроблоком. В этой схеме сформированы следующие критерии качества: T – время переходного процесса; y_max – максимальное значение выходной переменной; y'_max – максимальное значение производной выходной переменной. Критерий T сформирован в блоке Язык программирования, программа которого имеет вид

Листинг 7

input e;
if abs(e)>0.05 then T=time;
output T;

Рис. 9. а – модель оптимизируемой системы; б – объект (макроблок).

Потребуем, чтобы добротность по ускорению была не менее 4 с-2. Зададим ограничения на оптимизируемые параметры: 0 ≤ Kp ≤ 10, 4 ≤ Ki ≤ 20, 0 ≤ Kd ≤ 20. Зададим ограничения на критерии: 0 ≤ T ≤ 2, 0 ≤ y_max ≤ 1.3, 0 ≤ y'_max ≤ 3. Ограничения на параметры и критерии, а также метод оптимизации и вид (способ формирования) общего критерия (аддитивный, квадратичный, минимаксный, мультипликативный) задаются в диалоговом окне Параметры оптимизации.

Зададим в окне глобальных параметров начальные значения оптимизируемых параметров: Kp = 1; Ki = 5; Kd = 1. Выполнив моделирование, получим переходный процесс, показанный на рис. 10 красной линией. Такой процесс не удовлетворяет заданным требованиям, поэтому произведем оптимизацию, в результате которой получим Kp = 4.3; Ki = 4.0; Kd = 1.1. По желанию пользователя эти значения автоматически переносятся в окно глобальных параметров. Переходный процесс при таких параметрах (синяя линия на рис. 10) удовлетворяет всем заданным требованиям.

Рис. 10. Переходный процесс до оптимизации (красная линия) и после оптимизации (синяя линия).

4. Анализ

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

Расчет частотных характеристик выполняется путем подстановки s = jω и решения полученной системы линейных алгебраических уравнений с комплексными коэффициентами. При наличии в системе дискретных блоков принимается z = ejωT и используется гармонически линеаризованная модель экстраполятора нулевого порядка, задаваемая передаточной функцией W(s) = (1 - e-sT) / (sT), где T – период экстраполяции.

При расчете частотных характеристик необходимо задать с помощью блоков В память имена переменных, относительно которых будут рассчитаны характеристики. Далее эти имена используются для указания входов и выходов; кроме этого задается вид характеристики (АЧХ, ЛАХ, ФЧХ, Вещественная, Мнимая). Можно также построить годографы Найквиста, Попова, Михайлова, кривую D-разбиения. Результаты выводятся в соответствующее графическое окно.

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

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

5. Синтез

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

5.1. Частотный метод

Приведем основные соотношения частотного метода, реализованного в ПК «МВТУ». В общем случае структура линейной системы с одним входом и одним выходом имеет вид, показанный на рис. 11, где u0 – входное воздействие, u1 – сигнал управления, y0 – выходная (управляемая) переменная, y1, y2, ..., ym – измеряемые переменные, G(s) = [Gij(s)] – матричная передаточная функция неизменяемой части размером (m + 1) · 2, K(s) = [K1(s), K2(s), ..., Km(s),] – векторная передаточная функция регулятора. Примем передаточные функции регулятора в виде

(5.1)

Ki(s) = Ni(s) / D(s), i = 1, ..., m

где N1(s), ..., Nm(s), D(s) – полиномы числителей и общего знаменателя.

Рис. 11. Структура одномерной системы.

Обозначим через W(s) передаточную функцию системы от входа u0 к выходу y0, и пусть Wэт(s) – желаемая передаточная функция (эталонная модель), удовлетворяющая заданным требованиям к динамике системы. Зададим набор частот ω1, ω2, ..., ωn и потребуем, чтобы выполнялись условия

(5.2)

W(jωi) = Wэт(jωi), i = 1, ..., n.

На основе соотношений (5.2) можно составить систему линейных алгебраических уравнений, позволяющую найти неизвестные коэффициенты передаточных функций (5.1) регулятора. Хотя бы один из этих коэффициентов должен быть задан заранее, иначе система будет иметь бесконечное множество решений. Некоторые коэффициенты могут быть заданы из условий требуемого порядка астатизма, заданной добротности, фильтрующих свойств регулятора и т.д. Полученная система уравнений имеет комплексные коэффициенты, поэтому ее можно представить в виде системы из 2n уравнений с действительными коэффициентами. Пусть p – число неизвестных параметров регулятора, тогда из условия равенства числа уравнений и числа неизвестных получим при четном p число узловых частот n = p / 2. При нечетном p для узловой частоты ω1 формируется только одно уравнение по методу наименьших квадратов, тогда n = (p + 1) / 2.

При синтезе частотным методом в ПК "МВТУ" формируются две модели, одна из которой – модель синтезируемой системы, а вторая (эталонная) имеет желаемую передаточную функцию. В каждой из этих моделей с помощью блока В память обозначаются вход и выход. Обозначаются также входы и выход регулятора. В общем случае нужно найти передаточные функции регулятора (от каждого из его входов к выходу), в частном случае это могут быть коэффициенты. В специальном окне задаются: тип регулятора (непрерывный или дискретный), набор узловых частот, имена входов и выходов, порядки и некоторые из коэффициентов полиномов, входящих в (5.1).

Формирование желаемой передаточной функции по заданным требованиям является нетривиальной задачей, поэтому в ПК «МВТУ» создан каталог типовых эталонных моделей (\Demo\Синтез\Etalons). Необходимо только выбрать модель, соответствующую нужному типу ЛАХ и задать некоторые ее параметры (обычно это значения добротности и запаса устойчивости по фазе).

Для рассмотренной выше системы управления с ПИД-регулятором выберем эталонную модель с астатизмом 2-го порядка, в которой зададим добротность 5 с ‑2 и запас по фазе 45°. Узловые частоты зададим 0.1 и 3.5 с ‑1 (последнее значение примерно равно частоте среза эталонной системы). Подготовленная для решения задачи синтеза схема (файл \Demo\Синтез\Freq_Syn\Synt_PID.mrj) показана на рис. 12. Определяемыми параметрами являются здесь коэффициенты весового сумматора. В обратные связи включены ключи, которые используются для размыкания систем при расчете частотных характеристик. В результате синтеза получим Kp = 4.7; Ki = 5.0; Kd = 1.1. Частотные характеристики (ЛАХ и ФЧХ) эталонной и синтезированной систем приведены на рис. 13. Видно, что на низких и средних частотах кривые практически совпадают. Переходные характеристики в эталонной и синтезированной системах, представленные на рис. 14, также практически совпадают.

Рис. 12. Модель для решения задачи синтеза частотным методом.

Рис. 13. ЛАХ и ФЧХ эталонной системы (красная линия) и синтезированной системы (синяя линия)

Рис. 14. Переходная характеристика эталонной системы (красная линия) и синтезированной системы (синяя линия).

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

5.2. Корневой метод

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

Для рассматриваемой системы управления с ПИД-регулятором (файл \Demo\Синтез\ Root_Syn\Synt_PID.mrj) зададим полюсы -2 + j; - 2 - j; -3. В результате расчета получим Kp = 3.2; Ki = 2.4; Kd = 0.72 и полюсы -2 + j; - 2 - j; -3; -5.0; -35; -100. Заданные полюсы действительно оказались доминирующими. Если же задать полюсы ‑3; ‑4; ‑5, то один из полюсов синтезированной системы окажется равным ‑0.21, но тогда заданные полюсы уже не будут доминирующими. В таких случаях следует либо изменить структуру регулятора, либо задать другое расположение полюсов.

6. Исследование адаптивного ПИ-регулятора

ПК "МВТУ" является удобным инструментом не только для решения инженерных задач, но и при выполнении научных исследований. В частности, он использовался при разработке и отладке алгоритмов настройки и адаптации регуляторов на основе интерполяционного метода [14]. Рассмотрим применение ПК "МВТУ" для исследования адаптивной системы с ПИ‑регулятором.

Настройка параметров регулятора осуществляется на основе интерполяционного условия

(6.1)

,

где G(s) и Wэт(s) – передаточные функции объекта и разомкнутой эталонной системы, а частота настройки выбрана равной частоте среза эталонной системы ωс. Из (6.1) получим

(6.2)

.

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

Схема оценивания значений передаточных функций при комплексных значениях аргумента на основе взвешенного метода наименьших квадратов (МНК) была предложена в [14]. Эта схема использовалась при построении в ПК «МВТУ» моделей адаптивных систем с простейшими регуляторами (ПИ, ПД, ПИД). Проведенные эксперименты показали, что предложенная схема адаптации обеспечивает стабилизацию системы при изменении параметров объекта, но для настройки в заданной точке частотной характеристики входное воздействие должно иметь достаточно широкий спектр. В противном случае, например, при медленно изменяющемся входном воздействии, система стабилизируется, но не всегда настраивается на эталонную модель.

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

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

Рис. 15. а – модель адаптивной системы; б – регулятор; в – объект; г – устройство адаптации.

Модель адаптивной системы с ПИ-регулятором, сформированная в ПК "МВТУ" (файл \Demo\Адаптивные_регуляторы\Адап_ПИ.mrj) показана на рис. 15. На вход регулятора, кроме сигнала ошибки, поступает с выхода устройства адаптации векторный сигнал настраиваемых параметров Z = [Kp, Ki]. Объект представляет собой три последовательно соединенных апериодических звена с переменными коэффициентами. Модель объекта (рис. 14 в) состоит из двух блоков Язык программирования. В первом из них задается программа изменения параметров объекта. В нашем примере:

Листинг 8

if time<20 then K=1 else K=8;
if time<130 then T1=0.4 else T1=2;
T2=0.2; T3=0.1; P=[K,T1,T2,T3];
output P[4];

(time – модельное время). Вектор параметров P передается в блок, задающий уравнения объекта, а также (с помощью блока В память) в макроблок с подписью "Частотные хар-ки". Уравнения объекта задаются в виде

Листинг 9

input u,P[4];
init x1=0,x2=0,x3=0;
K=P[1]; T1=P[2]; T2=P[3]; T3=P[4];
x1'=(K*u-x1)/T1;
x2'=(x1-x2)/T2; x3'=(x2-x3)/T3;
output x3;

Устройство адаптации (рис. 7г) содержит в своем составе модель разомкнутой эталонной системы в виде передаточной функции sWэт(s). Примем

,

тогда параметры эталонной системы очень просто выражаются через задаваемые значения частоты среза ωс и запаса по фазе Δφ:

.

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

(6.3)

,

где β – скорость адаптации (подходящее значение – β = 0.1ωс). Передаточная функция H(s) = 1/s выбрана из условия H(jωс) = 1/(ωс), что обеспечивает выполнение соотношения (6.2) для оценки Kp. Отметим, что передаточные функции Wэт(s), F(s), H(s) могут быть и другими.

Оценивание параметров регулятора на основе соотношений (6.2) осуществляется в блоке "МНК" с помощью уравнений

Листинг 10

init Kp=Kp0,Ki=Ki0,v=v0;
Kp'= (u1-Kp*y1)*y1/v;
Ki'= (u2-Ki*y1)*y1/v;
v' = -beta*v + y1^2;
Z=[Kp,Ki];
output Z[2],v;

Здесь beta – значение β в (6.3). Вектор Z передается в регулятор, а также (с помощью блока В память) в макроблок "Частотные хар-ки", который обеспечивает расчет и отображение частотных характеристик, запаса по фазе и частоты среза. Переменная v используется только для отображения на график.

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

Листинг 11

dfi=60;                //запас по фазе, град
wc=2;                  //частота среза, рад/c
Kp0=0; Ki0=6; v0=1e-6; //начальные значения
beta=0.1*wc;           //скорость адаптации, 1/с
ctdfi=ctg(dfi*pi/180);
Kэт=wc*sqrt(1+ctdfi^2); Tэт=ctdfi/wc;

Задающее воздействие примем линейно нарастающим: g(t) = 0.01t. Такое воздействие – одно из наиболее неблагоприятных для адаптации. Процессы в адаптивной системе показаны на рис. 16. При заданных параметрах система неустойчива в начальный момент модельного времени, а также теряет устойчивость при t = 20, когда коэффициент K в модели объекта увеличивается в 8 раз. Система после этого быстро стабилизируется, но настройка на заданные значения запаса по фазе и частоты среза требует некоторого времени, в течение которого переменная v уменьшается до значения, при котором возникают устойчивые автоколебания. В установившемся режиме ошибка слежения не превышает 0.0056 (в эталонной системе установившаяся ошибка при таком воздействии равна 0.0043).

Рис. 16. Процессы в адаптивной системе.

Эксперименты, проведенные с рассмотренной моделью и аналогичными моделями, позволили сделать некоторые предположения относительно условий, при которых обеспечивается самонастройка адаптивной системы на заданные значения показателей качества. Многие из экспериментов проводились в режиме КОНТРОЛЬ И УПРАВЛЕНИЕ. Заданная настройка обеспечивалась во всех тех случаях, когда параметры объекта позволяли синтезировать систему заданного качества частотным методом. Таким образом, с помощью ПК "МВТУ" удалось показать перспективность предложенного способа построения адаптивных систем.

Заключение

Опыт решения большого числа тестовых и прикладных задач показал, что ПК «МВТУ» является удобным, эффективным и надежным средством для исследования и проектирования сложных технических систем, не уступающим по своим возможностям лучшим зарубежным аналогам. ПК «МВТУ» активно используется в учебном процессе МГТУ им. Н.Э. Баумана и ряда других университетов России и стран СНГ.

Программный комплекс внедрен в ряд проектных разработок предприятий Минатома России:

В рамках одной статьи невозможно показать все возможности ПК «МВТУ». Более подробную информацию, включая подробную техническую документацию, а также последние версии ПК «МВТУ» (демонстрационную и открытую учебную) с набором демонстрационных примеров, можно получить на сайте [1].

Литература и Интернет

  1. Сайт разработчиков ПК «МВТУ».
  2. Козлов О.С., Кондаков Д.Е., Скворцов Л.М. и др. Программный комплекс для исследования динамики и проектирования технических систем // Информационные технологии. 2005. № 9.
  3. Козлов О.С., Медведев В.С. Цифровое моделирование следящих приводов. В кн.: Следящие приводы: В 3‑х т. / Под ред. Б.К. Чемоданова. М.: Изд-во МГТУ им. Н.Э. Баумана, 1999. Т. 1. С. 711–806.
  4. Карташов Б.А., Карташов А.Б., Козлов О.С. и др. Практикум по автоматике. Математическое моделирование систем автоматического регулирования. М.: КолосС, 2004. 184 с.
  5. Марецкая В.В. Моделирование технологических процессов механической обработки с использованием программного комплекса «Моделирование в технических устройствах» («МВТУ») // Изв. вузов. Машиностроение. 2004. № 4. С. 39–52.
  6. Бенькович Е.С., Колесов Ю.Б., Сениченков Ю.Б. Практическое моделирование динамических систем. СПб.: БХВ-Петербург, 2002. 464 с.
  7. Кузин Д.А. Использование пакета Stateflow для создания модели силы трения // Exponenta Pro. Математика в приложениях. 2004. № 2. С. 10–16.
  8. Скворцов Л.М. Явные адаптивные методы численного решения жестких систем // Математическое моделирование. 2000. № 12. С. 97-107.
  9. Скворцов Л.М. Диагонально неявные FSAL-методы Рунге-Кутты для жестких и дифференциально-алгебраических систем // Математическое моделирование. 2002. Т. 14. № 2. С. 3–17.
  10. Козлов О. С., Скворцов Л. М. Тестовое сравнение решателей ОДУ системы MATLAB // Всероссийская научная конференция «Проектирование научных и инженерных приложений в среде MATLAB»: Труды. М.: Изд-во ИПУ РАН, 2002. С. 53–60.
  11. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.
  12. Крутько П. Д., Максимов А. И., Скворцов Л. М. Алгоритмы и программы проектирования автоматических систем. М.: Радио и связь, 1988. 306 с.
  13. Скворцов Л.М. Интерполяционные методы синтеза систем управления // Проблемы управления и информатики. 1998. № 6. С. 25–30.
  14. Скворцов Л.М. Интерполяционный метод автоматической настройки регуляторов // Изв. РАН. Теория и системы управления. 1998. № 6. С. 100–103.

15.06.2005