Моделирование систем управления

МОДЕЛИРОВАНИЕ СИСТЕМ УПРАВЛЕНИЯ 1. Два аспекта понятия моделирования. Понятие об идентификации. Моделирование в данном курсе понимается в двух аспектах: 1. Моделирование как процесс создания моделей. 2. Моделирование как процесс исследования свойств систем с использованием моделей. Идентификация  процедура получения математической модели по наблюдениям за входами и выходами объекта. Термин 'идентификация' возник от позднелатинского identifico, что означает равнозначность, соответствие. Для решения задачи управления всегда необходимо было изучить объект и получить его математическое описание. Создание такой модели и определение ее параметров и есть идентификация. В философском плане все наши представления об окружающем мире суть модели. Каждая модель 'хромает'. Всякая модель отражает объект с какой-либо стороны. Совершенствовать модель можно до бесконечности. Критерий истинности  всегда практика. Формирование моделей на основе результатов наблюдений и исследование их свойств  вот, по существу, основное содержание науки. Решение задачи построения математических моделей динамических систем по данным наблюдений за их поведением составляет предмет теории идентификации, которая тем самым становится элементом общей научной методологии. А поскольку мы 'окружены' динамическими системами, методы идентификации систем имеют широкие приложения. 2. Причины необходимости создания новых моделей Теория автоматического управления прошла путь от систем стабилизации (регулятор Уатта) до систем управления ТЭП (техникоэкономическими показателями) процессов. Развилась также и техника управления. С развитием вычислительной техники появилась возможность управления сложными технологическими процессами на новом качественном уровне (с применением микропроцессорных устройств), что требует:  более тщательного изучения объекта управления (ОУ),  исследование более сложных ОУ (увеличение их числа),  исследование ОУ в более широком диапазоне изменения переменных и т.д. Все это обуславливает развитие методов идентификации и разработку теории идентификации. 2 На современном этапе требуется идентифицировать как устройства, так и технологические процессы. С точки зрения теории и практики идентификации устройства и процессы представляют собой динамические системы со входами и выходами, но процессы представляют больший интерес. 3. Характеристики объектов и процессов, которые надо учитывать при создании моделей Рассмотрим основные особенности непрерывных технологических процессов, которые необходимо учитывать при построении математических моделей объектов управления: 1. Наличие большого числа входных и выходных переменных. Выходные переменные представляют собой параметры, входящие в критерий оптимальности, или определяющие качество переходных процессов. Входные воздействия можно разделить на управляющие и возмущающие. Последние могут быть контролируемыми или неконтролируемыми, причем точное число неконтролируемых переменных неизвестно. 2. Наличие физических и технологических развязок. Физические развязки связаны с тем, что воздействия одной физической природы не влияют на параметры другой физической природы. Технологические развязки представляют собой различные бункеры, накопительные емкости, сборники и т.п. Благодаря технологическим и физическим развязкам, процесс может быть представлен как совокупность отдельных объектов управления. 3. Существенность динамических свойств объектов управления. Значения выходных переменных в данный момент времени зависит не только от настоящих, но и от предыдущих входных воздействий. Время памяти является конечным и для различных технологических процессов может измеряться от нескольких секунд до десятков часов. Изменения выходных параметров под воздействием управляющих переменных должны протекать быстрее, чем под воздействием возмущений. 4. Нелинейность объектов управления. Как правило, для более производительного протекания процессов необходимо варьировать управления в более широких пределах, что приводит к проявлению нелинейностей присущих всегда реальным процессам. При малых изменениях переменных от номинального значения объекты можно линеаризовать. 5. Нестационарность статических и динамических характеристик процессов. 3 Обычно нестационарность связана с изменением условий работы оборудования, изменением режимов технологических процессов, действием неучтенных возмущений. 6. Неконтролируемость ряда выходных переменных технологических процессов из-за отсутствия необходимых измерительных преобразователей. Существуют также процессы с редко контролируемыми параметрами, например, по результатам лабораторных анализов, проводимых периодически. 7. Высокий уровень помех при измерении входных и выходных переменных, связанный с несовершенством существующих измерительных преобразователей. Высокий уровень помех обусловлен в значительной степени интенсификацией производства. 4. Приемы упрощения моделей При построении математического описания технологических процессов стремятся получить наиболее простые модели. При этом пользуются следующими приемами упрощения: 1) расчленение сложной системы на ряд более простых, т.е. декомпозиция; 2) выделение существенных связей, учет несущественных воздействий в виде переменных параметров модели; 3) пренебрежение динамическими свойствами и получение на первом этапе статических моделей; 4) линеаризация нелинейных процессов в определенной области изменения переменных; 5) приведение процессов c распределенными параметрами к процессам с сосредоточенными параметрами; 6) замена нестационарных процессов кусочно-стационарными. 5. Этапы построения моделей Выделим пять этапов построения моделей: 1) определение цели получения модели; 2) определение ограничений и условий, учитываемых построении моделей; 3) выбор подхода к решению задачи получения модели; 4) определение класса модели; 5) выбор метода 6. при Определение цели получения модели Один и тот же физический объект можно рассматривать с точки зрения решения различных задач. Вначале надо определить, для какой цели требуется модель объекта. Цели получения модели могут быть, например, следующие: 4 1) 1) 2) 3) 4) 5) расчет не измеряемых параметров; контроль состояния объекта; прогнозирование параметров объекта; управление процессом в статике или динамике; исследование системы управления объектом; анализ и исследование отдельных особенностей объекта управления и др. Достижение ряда целей не требует полного описания объекта, достаточно иметь описание лишь тех частных качеств, которые требуются для достижения поставленных целей. Гомоморфная модель  модель, отражающая отдельные частные свойства объекта или процесса. Изоморфная модель  полная модель объекта. Исчерпывающее описание объекта в виде изоморфной модели возможно при наличии полной информации об объекте. Последнее, во-первых, дорого, а во-вторых, достижимо, как правило, только теоретически. Модели одного и того же физического процесса могут иметь мало общего, если они разрабатывались для разных целей. 7. Определение ограничений и условий, учитываемых при построении моделей Ограничения, учитываемые при построении модели делятся на: 1) технико-экономические и 2) технические. Технико-экономические ограничения определяются техническими требованиями, экономическими показателями, заданными при проектировании системы управления. К их числу относятся: 1) затраты на получение модели (финансовые и временные); 2) требуемая точность модели; 3) возможность уточнения и исправления модели в процессе эксплуатации системы управления; 4) техническая обеспеченность системы управления (наличие необходимых измерительных преобразователей, мощность вычислительного комплекса и т.п.). Технические условия определены природой изучаемого объекта управления. Во-первых, это условия, указывающие на характер исследуемого объекта: 1) линейность, 2) стационарность, 3) инерцыонность, 4) размерность и т.д.; во-вторых, условия, зависящие от вида входных и выходных сигналов: 1) эргодичность, 5 2) 3) 4) 5) наконец, объекта: 1) 2) 3) 4) стационарность, уровень помех, место приложения помех, характер помех; можно выделить условия, которые определяются режимом работы динамический, статический, допускающий активный эксперимент, не допускающий активный эксперимент. 8. Выбор подхода к решению задачи получения модели Далее необходимо выбрать подход к получению математической модели: 1) аналитический, 2) экспериментальный, 3) экспериментально-аналитический. Аналитический подход заключается в аналитическом выводе уравнения объекта, исходя из основных физических и химических законов, и определение входящих в них параметров по опытным данным, полученным из специально поставленных экспериментов. Этот подход называют иногда моделированием. Экспериментальный подход имеет ряд названий: кибернетический, метод черного ящика, метод идентификации. Идентификацией называют процесс построения математической модели технического объекта по его измеряемым входным и выходным сигналам. При решении задачи идентификации в широком смысле начальная информация об объекте отсутствует и система представляется в виде "черного ящика". Если структура объекта и, следовательно, его модели известна, а по результатам эксперимента определяются параметры модели, то такой подход называется экспериментально-аналитическим или идентификацией в узком смысле. Этот подход называют также методом "серого ящика" или параметрической идентификацией. Аналитический подход 1) требует обычно больших затрат времени, 2) дает модель сравнительно невысокой точности, 3) дает модель работоспособную в широком интервале изменения переменных. Экспериментальный подход 1) позволяет быстро находить модели, 2) дает модель, довольно точно описывающую отдельные режимы работы объекта, но 6 3) работоспособность модели ограничена исследованного изменения переменных. 9. узкими интервалами Определение класса модели. Выбор метода решения задачи и ее решение Класс модели определяется математической записью модели 1) параметрической или 2) непараметрической. Например, дифференциальное уравнение  это параметрическая форма, а модель в виде импульсной характеристики  непараметрическая форма. Выбор метода решения задачи идентификации связан со знанием классификации их. Методы идентификации можно классифицировать по ряду признаков: а) по характеру объекта управления это методы идентификации: 1) в статике или динамике; 2) линейных или нелинейных ОУ; 3) стационарных или нестационарных; 4) одномерных или многомерных ОУ; 5) ОУ с сосредоточенными или с распределенными параметрами; 6) непрерывных и дискретных; 7) детерминированных или стохастических. Наиболее детально разработаны методы идентификации статических и динамических линейных стационарных объектов с сосредоточенными параметрами. б) по характеру входного воздействия различают 1) методы активного эксперимента и 2) методы пассивного эксперимента. в) по степени оперативности различают такие методы как: 1) оперативная идентификация (в темпе поступающей информации) и 2) ретроспективная идентификация (послеопытная). Большинство методов идентификации являются ретроспективными. Оперативные методы используют рекуррентные алгоритмы, позволяющие уточнять модель на каждом шаге работы алгоритма. Общая схема процесса оперативной идентификации объекта управления показана на рис. 1. 7 x  помеха на входе y Объект управления  x' Модель объекта y'   помеха на выходе Формирование критерия близости коррекция модели Минимизация критерия Рис.1. Схема оперативной идентификации Возможно различать методы также в зависимости от вида критериев близости модели и объекта. С математической точки зрения задача идентификации представляет собой задачу на поиск экстремума, условного или безусловного, функции многих переменных для параметрической идентификации в статике или на экстремум функционала для непараметрической идентификации или для задач динамики. С другой стороны задача идентификации относится к обратным задачам, которые, как правило, труднее прямых задач и часто приводят к т.н. некорректным задачам. Таким образом, в результате решения задачи идентификации строится математическая модель в заданном классе, выбранной формы, оптимальная в смысле некоторого критерия близости при исходной информации в виде результатов измерения входных и выходных сигналов объекта управления. При разработке моделей физического процесса наряду с классическими методами наблюдения, систематизации, построения гипотез и испытаний видное место занимают интуиция и здравый смысл. Интуиция играет важную роль при формировании основных допущений в установлении основных зависимостей между ключевыми переменными, а также при выработке первоначального подхода к построению модели физического процесса. Здравый смысл требуется для обеспечения равновесия между точностью и полнотой описания модели, с одной стороны, и сложностью и высокой стоимостью управления, с другой. 8 10. Общие принципы построения аналитических моделей Для получения математического описания объектов управления не всегда возможно проведение экспериментального исследования (например, при проектировании новых производств). При отсутствии возможности экспериментальных исследований математическое описание объекта управления может быть получено с помощью аналитических методов (если для этого есть условия, т.е. в первую очередь, известно с теоретической точки зрения содержание процессов, происходящих в объекте управления). В некоторых случаях (например, при составлении моделей механических систем роботов) теория процессов, происходящих в объекте известна настолько хорошо, что проще аналитически вывести математическую модель, чем предпринимать дорогостоящее экспериментальное исследование. Для составления математического описания стараются разделить исследуемый объект на более простые составляющие (звенья), используя, в частности, различные физические или технологические развязки. Затем по математическим описаниям отдельных звеньев определяется модель объекта в целом. Общего метода, пригодного во всех случаях, для аналитического вывода модели не существует. В каждой области техники используют различные приемы, исторически сложившиеся при моделировании объектов из данной области. Разумеется, что общим для всех приемов является использование единых физических законов. Наиболее часто для составления математического описания объектов управления используются уравнения материального и энергетического баланса. Уравнение материального баланса основано на законе сохранения массы вещества и может быть записано в следующей форме: n  G i  0, i 1 где Gi – весовой расход [кг/сек] i-го вещества. Расходы веществ, поступающих в объект, считаются положительными, а выходящих из него – отрицательными. Записанное уравнение справедливо для установившегося состояния, т.е. оно является уравнением статики. Для описания динамики подобное уравнение составляется для бесконечно малого промежутка времени и добавляется бесконечно малое изменение массы, приводящее к несбалансированности масс в динамике (которая, в конечном счете, и является причиной движения динамической системы) к балансу. Эта добавка позволяет, таким образом, описать динамику в виде уравнения 9 n  dG i  dG. i 1 Часто для получения математического описания объектов управления используются уравнения энергетического баланса. В установившемся режиме количество энергии, переходящей в объект, равно количеству энергии, уходящей из него. Если подводимую энергию считать положительной, а уходящую – отрицательной, то энергетический баланс для установившегося режима запишется в виде n  Wi  0, i 1 где Wi – i-й поток энергии. Динамическое уравнение теплового баланса имеет вид: приход – расход = накопление или dWвых Wвх  Wвых  k , dt где Wвх – входной поток энергии, Wвых – выходной поток энергии, k – коэффициент. Анализ разнообразных процессов, протекающих в химикотехнологических аппаратах, позволяет выделить некоторую совокупность типовых объектов управления, описываемых с помощью типовых математических моделей. Перечень и характеристики типовых моделей приводятся в специальной литературе [Кафаров, Методы кибернетики, 1971]. При аналитическом составлении моделей механических систем часто используются уравнения баланса сил в статике, а в динамике в эти уравнения вводятся силы инерции (принцип Даламбера). В электротехнике, гидравлике, теплотехнике, газовой динамике, пневматике, химической кинетике и др. установились свои приемы описания динамики объектов. Каждый из этих приемов требует изучения, хотя все они вытекают из фундаментальных законов сохранения вещества и энергии. Рассмотрим применение аналитического подхода на нескольких конкретных примерах. 11. Модель поплавкового уровнемера Пусть из исследуемого объекта выделено звено, функцией которого является преобразование уровня жидкости в перемещение пластины емкостного преобразователя. 10 y x=h V - объем поплавка +  Рис. 2. Поплавковый преобразователь уровня Как видно из рисунка, звено состоит из сферического поплавка, погруженного в исходном положении наполовину в жидкость и с помощью штока прикрепленного к левому плечу рычага. Конец правого плеча рычага является пластиной конденсатора емкостного преобразователя. Известно, что электростатическими силами, действующими со стороны емкостного преобразователя на динамику перемещения поплавка можно пренебречь ввиду их малости. Это пример физической развязки, которая позволяет выделить поплавок с рычагом в отдельное звено. Входным сигналом звена будем считать изменение уровня жидкости h относительно некоторого исходного, условно принятого за 0. Положительным будем считать значения h при увеличении уровня жидкости. За выходной сигнал для удобства примем не перемещение пластины емкостного преобразователя, а перемещение, обозначенное как y, штока поплавка. При равенстве плеч рычага и при абсолютной жесткости рычага эти перемещения будут равны. Таким образом, анализируемое звено можно представить в виде блока с одним входом и выходом x=h W(p) y Рис. 3. Искомая модель поплавкового преобразователя Исследуем свободное движение поплавка. Для этого сообщим системе некоторое начальное смещение h0 = y0 и начальную скорость y (0) и, убрав в некоторый момент t = 0 внешние силы (которые нужны были для задания начальных значений), будем наблюдать свободное движение (в данном случае механическое движение) системы. В любой момент времени на систему действуют: выталкивающая сила Архимеда Pa(y), направленная вверх (примем это направление за положительное); сила веса, направленная вниз и равная mg (m  масса 11 системы, g  ускорение свободного падения); сила трения, направленная против движения и пропорциональная (как следует из изучения физики движения тел в жидкости при малых скоростях) скорости движения, т.е. равная kdy/dt; сила инерции, которую необходимо учесть в динамике и которая по принципу Даламбера направлена против движения и равна md2y/dt2. Таким образом, уравнение баланса сил имеет вид: dy d2y m  0. (*) dt dt 2 В установившемся состоянии y = 0, dy/dt = 0, d2y/dt2 = 0 и из уравнения динамики получим уравнение статики Pa ( y)  mg  k Pa(0) = mg. (**) Рассмотрим выталкивающую силу и ее зависимость от y. По закону Архимеда можно записать, что она равна объему погруженной части поплавка Vп(y), умноженной на плотность жидкости d и на ускорение свободного падения g, т.е. равна весу вытесненной жидкости: Pa(y) = V(y)dg. Обозначим полный объем поплавка радиуса R через V = 4R3/3 и будем считать, что в установившемся состоянии поплавок погружен наполовину и при смещении от этого положения изменения погруженной части происходят за счет малых (при малых y) изменений объема средней части поплавка. Будем считать, что в средней части поплавок имеет приближенно цилиндрическую форму с площадью основания R2. Тогда Vп(y) = V/2  R2y; Pa(y) = (V/2  R2y)dg Т.к. (V/2)dg = Pa(0), то Pa(y) = Pa(0)  R2dgy. Подставляя в (*), получаем: d2y dy  R 2dgy  0. dt dt Разделив обе части уравнения на коэффициент при y, получим математическую модель в хорошо известной форме: m T2 2 k d2y dt 2  2T dy  y  0, dt 12 где T  m 2 R dg  постоянная времени звена;   k 2 mR 2dg  степень затухания. Преобразовав последнее уравнение по Лапласу, можно представить искомую модель звена в виде: x 1 T2p2 + 2Tp + 1 y Рис. 4. Модель поплавкового преобразователя Рассмотрим параметры звена T и . V 4 1 Учитывая, что Pa (0)  mg  dg  R 3 dg , получаем: 2 3 2 4R 3 d. 6 Подставив это выражение для массы в выражение для T, получим: 2R T . 3g m Последнее выражение позволяет вычислить постоянную времени (как важнейшую характеристику динамики), если известен конструктивный параметр R (например, при R = 1 м и g = 10 м/с T = 0.26 с), или вычислить конструктивный параметр R, задавшись значением T. Что касается второго параметра  степени успокоения , то он зависит от коэффициента трения k, численное значение которого теоретически получить не удается. Для определения этого параметра необходимо прибегнуть к специально поставленному эксперименту. Необходимо практически изготовить опытный поплавок из материала, который будет использоваться при изготовлении реального устройства, и записать свободные затухающие колебания поплавка в той же жидкости, которая будет в реальном объекте. Используя затем методы идентификации (изучаемые ниже), по полученной экспериментальной записи можно получить численные значения  или k. Этот процесс аналитического вывода модели и экспериментального подхода для определения параметров является характерным для большинства случаев т.н. аналитического подхода для получения математического описания объекта. В связи с этим часто утверждают, что нет чисто аналитического подхода, а есть аналитическо-экспериментальный подход. 13 Отметим (на данном характерном примере), что аналитический вывод позволяет получить пригодные результаты при условии введения определенного числа упрощений и допущений. Например, в данном случае не учтены такие эффекты как:  трение о воздух и связанные с этим потери энергии,  влияние температуры на линейное расширение материалов,  влияние соотношения объемов поплавка и емкости с жидкостью,  не идеальность поверхности жидкости и ряд других. При учете всех факторов модель может получиться слишком сложной и в то же время малопригодной для дальнейшего решения задачи управления. 12. Модель процесса теплопередачи Рассмотрим аналитическое описание процессов передачи тепла. Передача тепла относится к довольно сложным задачам описания процессов в форме дифференциальных уравнений. Это связано в первую очередь с тем, что передача тепла осуществляется на основе трех различных физических процессов: теплопередачи, когда на молекулярном уровне быстродвижущиеся молекулы в точке с более высокой температурой за счет соударения с соседними молекулами отдают им часть своей энергии, увеличивая тем самым скорость медленных молекул и повышая в конечном счете температуру в соседней области; излучения и конвекции. Принятие во внимание сразу всех трех физических эффектов для описания передачи тепла обычно не используется. В каждом конкретном случае с самого начала принимается решение о преобладании того или иного эффекта и два других не рассматриваются. Рассмотрим процесс теплопередачи. Обычно задача ставится следующим образом: найти уравнения, связывающие изменение во времени температуры тела, находящегося в контакте с горячей средой, при заданном изменении температуры среды. В общем случае это приводит к дифференциальным уравнениям в частных производных, так как температура тела зависит не только от времени, но и от трех пространственных координат. Для упрощения задачи и сведения уравнений к обыкновенным дифференциальным уравнениям принимают ряд упрощающих допущений таких как одинаковость температуры горячей среды во всех точках, постоянство тепловых параметров и т.д. Кроме того, сводят задачу или к плоскому распространению тепла без учета краевых эффектов или используют какой-либо тип симметрии, наиболее подходящий к конкретной 14 ситуации, и тогда рассматривают распространение тепла вдоль оси симметрии. Рассмотрим для примера плоскую задачу передачи тепла от горячей среды (топки) с температурой Тс (температура среды) твердому телу, температуру которого обозначим Тn (температура приемника). Требуется построить модель x = Tc y = Tп W(p) Рис. 5. Искомая модель теплопередачи Будем рассматривать процесс в трубке с поперечным сечением S Тс Тп Рис. 6. Поток тепла от среды к приемнику Из физики известно, что количество тепла Q1, переданного от горячей среды холодному телу за время dt, пропорционально dt, площади S и разности температур Тс  Тn Q1 = c S (Тс  Тn) dt, где c  коэффициент теплопередачи, характеризующий данное тело (размерность его кал/(м2 град сек)). Количество тепла Q2, полученное приемником, пропорционально массе тела m и изменению температуры приемника dTn за то же время Q2 = C m dTn, где C  удельная теплоемкость (кал/(кг град) ). Приравнивая эти два количества тепла, получаем или: С m dTn = c S (Tc  Тn) dt C m dTn  Tn  Tc или Ty  y  x ,  cS dt где T  C m  постоянная времени модели в виде звена первого порядка:  cS 15 x = Tc 1 Tp + 1 y = Tп Рис. 7. Модель теплопередачи Полученная упрощенная модель позволяет в целом верно оценить инерционность звена, как объекта управления, и анализировать зависимость этой инерционности от массы, размеров и теплофизических свойств тела. Если числовые значения коэффициентов C и c отсутствуют в справочнике, то возможно потребуются специально поставленные эксперименты для их определения. Рассмотрим дополнительно ситуацию, когда два тела с разными теплофизическими свойствами находятся в контакте друг с другом и с горячей средой (теплообменник или датчик температуры в защитном чехле и т.п.). Тс Тп1 Тп2 Рис. 8. Теплопередача двум приемникам В этом случае задачу приближенно рассматривают как передачу тепла от горячей среды первому телу, которое непосредственно в контакте со средой, а затем рассматривают передачу тепла от первого тела ко второму. Если при малой толщине первого тела можно принять гипотезу, что его температура примерно одинакова в любой точке, то вторая задача решается аналогично первой и тогда получаем: T1y 1  y1  x , T2 y 2  y 2  y1, где y1  температура первого тела; y2  температура второго тела; T1, T2  постоянные времени, обусловленные первым и вторым телом соответственно. Исключая из этих уравнений переменную y1, получим модель в виде: T1T2 y 2  (T1  T2 ) y 2  y 2  x. Легко можно показать, что полученное дифференциальное уравнение второго порядка не может иметь колебательного решения. Это 16 обстоятельство соответствует реально наблюдаемым процессам теплопередачи и может служить в некоторой мере основанием допустимости принятых упрощений. Если бы конечное уравнение допускало колебательные решения, которые нельзя получить на практике, то скорей всего мы бы стали искать возможность построить более точную модель. 13. Модель смесителя. В смеситель (рис. 9) поступает два потока веществ А и B: поток QA [м3/сек] с концентрацией СA [кг/м3] и поток QB [м3/сек] с концентрацией СB [кг/м3]. На выходе смесителя поток QC c концентрацией смеси СC. мотор мешалки A B QA CA QB CB объем смесителя VC = const QC CC C Рис. 9. Смеситель двух компонентов Требуется построить модель x1 = CA x2 = CB W(p) y = CC Рис. 10. Искомая модель смесителя Примем допущения, что перемешивание мгновенное по всему объему и объем смесителя VC постоянен: VC = const. Последнее условие эквивалентно условию: QC = QA + QB. Материальный баланс в переходном процессе может быть выражен следующим образом: изменение количества вещества, находящегося в 17 растворе емкости, равно разности между количеством притекающего и вытекающего из емкости вещества за одно и то же время: dmC = QACAdt + QBCBdt  QCCCdt. (*) Масса вещества смеси mC = VC CC При постоянстве объема смеси изменение массы может быть только за счет изменения концентрации, т.е. dmC = VCdCC. Из (*) и (**) получаем dCC VC  Q C CC  Q A C A  Q BC B . dt Откуда имеем Ty  y  K1x1  K 2 x 2 , (**) где T Vc Vc QA ; ; K1   Qc Q A  QB QA  QB K2  QB . QA  QB Искомая модель может быть представлена, таким образом, в виде: x1 K1 1 Tp + 1 x2 y K2 Рис. 11. Структура модели смесителя Во многих случаях такой упрощенной модели вполне достаточно для управления процессом смешения. Более сложные модели получаются при учете непостоянства объема смеси, неодинаковости концентрации смеси по объему смесителя, при учете явлений кавитации и т.д. Усложнение модели может привести к тому, что она станет непригодной для решения задачи управления. 14. Модель реактора В отраслях химической промышленности распространены разнообразные реакторы с перемешивающими устройствами. Они 18 отличаются друг от друга по конструкции и по типу проходящих в них реакций. На рис. 12 показан реактор с мешалкой, в котором происходит реакция, в результате которой вещество А превращается в вещество В. мотор мешалки A Q CA объем реактора Vр = const Q CрА СрВ Рис. 12. Химический реактор Составим математическое описание объекта управления в предположении, что объем раствора VP [m3] в реакторе постоянен, т.е. расходы на входе и выходе равны (Q [m3/s]). x = CA [kg/m3] W(p) y1 = CрА [kg/m3] y2 = CрB [kg/m3] Рис. 13. Искомая модель реактора В искомой модели предусмотрено два выхода на тот случай, что управление может осуществляться или концентрацией вещества А в растворе (СPA), если цель состоит в уничтожении вещества А (переводом его в вещество В), или концентрацией вещества В, если получение его является целью реакции. Для составления математического описания необходимо в первую очередь знать динамику самого превращения вещества А в вещество В. Из химических исследований известно, например, что многие подобные реакции могут быть описаны постоянной величиной, называемой скоростью реакции. Скорость реакции k это количество вновь возникающего вещества в единицу времени, отнесенное к количеству вступившего в реакцию исходного вещества m / mA k B . (*) dt Возможны реакции, которые подчиняются другим закономерностям превращения. 19 Составим уравнения материального баланса отдельно для вещества А и вещества В. Для вещества А: Для вещества В: dmA = Q CA dt  Q CPA dt  mB. (**) dmB = mB  Q CPB dt. (***) Массу mB вещества В, возникшего за время dt, запишем как mB = k mA dt Кроме того, при условии идеального перемешивания можно записать и, следовательно: mA = VP CPA; mB = VP CPB dmA = VP dCPA; dmB = VP dCPB. Подставим последние соотношения в уравнения (**) и (***) VP dCPA = Q CA dt  Q CPA dt  k VP CPA dt VP dCPB = k VP CPA dt  Q CPB dt Откуда находим dy VP Q T1 1  y1  k1x; T1  ; k1  ; dt Q  kVP Q  kVP dy 2 V kV  y 2  k 2 y1; T2  P ; k 2  P ; dt Q Q В соответствии с полученными уравнениями структура модели реактора будет иметь вид, показанный на рис.14. T2 y1 x 1 T1p + 1 1 T2p + 1 y2 Рис. 14. Структура модели реактора Исключая из двух последних переменную y1, можно получить дифференциальных уравнений 20 d 2 y2 dy 2  y 2  k1k 2 x. dt dt С полученной моделью и коэффициентами можно выполнить массу различных анализов и преобразований, но постоянно надо иметь ввиду те допущения и упрощения, которые были приняты при выводе модели. Некоторые выводы могут прийти в противоречие с практикой. T1T2 2  (T1  T2 ) 15. Модель емкости с изменяющимся уровнем Пусть требуется построить модель для емкости, в которую подается жидкость и из которой вытекает жидкость. Входной переменной пусть будет расход на входе, а выходной уровень жидкости в емкости. Q1 [m3/s] p1 [N/m2] h V [m3] S [m2] Q2 p2 [N/m2] Рис. 15. Емкость с меняющимся уровнем Материальный баланс в данном случае может быть сведен к следующему: изменение объема жидкости за время dt равно разности объемов втекающей и вытекающей жидкости, т.е. или dV = Q1 dt  Q2 dt dV/dt = Q1  Q2. dV dh S . dt dt Истечение жидкости из емкости описывается уравнением Бернулли: Объем жидкости в емкости V = S h; откуда v 2  v 02 p p  h  h0  1 2 ,  2g где v  скорость истечения из сливного отверстия, v0  скорость изменения уровня в резервуаре, 21 h  h0  перепад высот, p1, p2  давление над жидкостью и за сливным отверстием,   удельный вес. При v >> v0; h0 = 0 и p1 = p2 имеем: v2  h; 2g v  2gh . Пусть S0  сечение сливного отверстия, тогда в идеальном случае Q2 = v S0. Реально последнее выражение записывается с некоторым коэффициентом , отражающим зависимость расхода Q2 от неравномерности скорости v, от качества обработки отверстия, от его формы и т.д. Q 2  S0 2gh . В результате получаем модель dh S  S0 2gh  Q1. dt Пример наглядно показывает, что даже в простой ситуации модель получается в виде нелинейного дифференциального уравнения. 16. Метод наименьших квадратов в одномерном случае Метод наименьших квадратов известен давно и применяется весьма широко и повсеместно. Преимущество метода  простые вычисления для определения степенной аппроксимации. Рассмотрим простейший случай. x Объект y Рис. 16. Одномерная система Проводим эксперимент N раз и получаем пары чисел xi, yi; i = 1, 2, ..., N 22 y y y0 x0 x x Рис. 17. Экспериментальные данные 1) Делается допущение, что y зависит от x линейно, а несоответствие линейной зависимости обуславливается случайными факторами. Т.о. ставится произвольно задача найти линейную зависимость ŷ i  a 0  a1x i с постоянными коэффициентами a0, a1, которые позволят для любого экспериментального значения xi вычислить некоторую оценку ŷ i , которая нас удовлетворит в каком-либо смысле. 2) Второе допущение состоит в том, что положительные и отрицательные отклонения должны компенсировать друг друга и искомая зависимость должна проходить через точку с координатами, равными средним значениям 1 N 1 N x   x i ; y   y i, N i 1 N i 1 т.е. искомое уравнение должно удовлетворять соотношению y  a 0  a1x. (*) 3) Следующее допущение касается аддитивности ошибки, т.е. мы принимаем, что y i  ŷ i   i  a 0  a1x i   i . (**) Перенесем начало координат в точку (центрированные) значения переменных y 0i  yi  y; ( x ,y ) x 0i  x i  x . и введем новые 23 Вычитая теперь из (**) уравнение (*), получим y 0i  a1x 0i  i ;  i  y 0i  a1x 0i . Обратим внимание на то, что все отмеченные выше допущения принимаются как аксиомы и направлены на то, чтобы получить ошибку как линейную функцию неизвестных коэффициентов аппроксимирующей зависимости i = f(a1). 4) Наконец, самое основное допущение оптимизации, который принимается в виде касается критерия N J   i2  min . i 1 Принятие критерия в таком виде (который дает имя методу) позволяет однозначно и легко вычислить искомый коэффициент известным методом поиска экстремума функции одной переменной: dJ d d   ( y0i  a1x 0i )2   ( y0i  a1x 0i )2  da1 da1 da1   2( y0i  a1x 0i )( x 0i )  0;  0i 0i . 2  y 0i x 0i  a1  x 0i  0  a1  2  x 0i y x Далее из (*) находится а0: a 0  y  a1x. 17. Метод наименьших квадратов в многомерном случае Будем рассматривать здесь многомерную модель как модель со многими входами и одним выходом: x1(t) x2(t) Система y(t) xk(t) Рис. 18. Многомерная система Сигналы на входах и выходах считаем центрированными, т.е. M[xj(t)] = 0, j = 1, 2, ..., k; M[y(t)] = 0. Модель ищем в виде линейной комбинации входных сигналов 24 ŷ( t )  b1x1( t )  b 2 x 2 ( t )  b k x k ( t )  k  b jx j ( t ). j1 Помеху считаем аддитивной y( t )  ŷ( t )  ( t )  k  b jx j ( t )  ( t ). j1 Выбираем N >> k моментов времени ti, i = 1, 2, ..., N и в каждый из этих моментов наблюдаем все входы и выход. Интервал между ti не играет роли. Важно только, что в момент взятия отсчетов система должна быть в установившемся состоянии, т.е. интервал должен быть достаточным для практического окончания переходных процессов при изменении очередного состояния системы. Введем обозначения: xj(ti) = xij; y(ti) = yi; (ti) = i. тогда y i   b jx ij   i ;i  1, N. j В матричной форме Y = XB + E, где y1 y Y 2 ;  yN x11 x X  21  x N1 x12 x 22  x N2     x1k x 2k ; B   x Nk b1 b2 ; E   bk 1 2 .  N Критерий приближения J   i2 i 1   1  2   N 2  E T E  N Так как E = Y - X B, то J = (Y - XB)T(Y - XB) = (YT – BTXT)(Y – XB) = = YTY – BTXTY – YTXB + BTXTXB Поскольку все слагаемые в последнем уравнении скаляры, то BTXTY = (BTXTY)T = YTXB Тогда J = YTY – 2BTXTY + BTXTXB 25 Вектор коэффициентов B можно вычислить из условия: dJ/dB = 0: dJ d(Y T Y  2B T XT Y  B T XT X  B)   dB dB d(B T XT Y) d(B T X T X  B) d( B T X T X  B )   2 X T Y  dB dB dB Здесь использовано правило для двух векторов y и z:  2 d(z T y ) y dz Далее воспользуемся правилом: если Q – квадратная матрица, а z – вектор соответствующего размера, то d(z T Qz )  (Q  Q T )z dz Для квадратной матрицы XTX и вектора B получаем:   d(B T XT X  B)  X T X  ( X T X) T B  2 X T X  B dB Тогда dJ/dB = -2XTY + 2XTXB = 0 Откуда получаем т.н. нормальные уравнения МНК: XTXB = XTY. Решение нормальных уравнений B = (XTX)-1 XTY Обычно приведенное решение используется в теоретическом анализе, а вычисление вектора коэффициентов B выполняется решением системы нормальных уравнений (например, методом Гаусса). Запишем матрицы системы x11 x 21  x N1 x11 x12  x1k x T X X  12 x 22  x N 2 x 21 x 22  x 2 k          x1k x 2k  x Nk x N1 x N 2  x Nk 26 2  x i1  x i1x i 2   x i1x ik 2   x i 2 x i1  x i 2   x i 2 x ik ; dimX T X  k  k;     2  x x x x  ik i1  ik i 2  x ik x11 x X Y  12  x1k T x 21 x 22  x 2k     x N1 x N2  x Nk y1  x i1yi y 2   x i 2 yi ; dimX T Y  k  1.  y N  x ik yi В одномерном случае (k = 1) положим xi1 = xi, тогда yi = b1 xi + i X T X   x i2 ; ( X T X) 1  1  x i2 b1  ( X T X) 1 X T Y  ; XT Y   x i yi ;  x i yi ,  x i2 что совпадает с ранее полученным ранее выражением. 18. Рекуррентный метод наименьших квадратов Пусть по N наблюдениям вычислен вектор коэффициентов BN = (XTX)-1 XTY = PN XTY где PN = (XTX)-1, и пусть в момент tN+1 получена еще одна серия отсчетов yN+1, xN+1,1, xN+1,2, ..., xN+1,k. Необходимо найти способ откорректировать вектор BN по информации, полученной на (N+1)-ом шаге, т.е. найти способ вычисления по формуле BN+1 = BN + f(yN+1, xN+1,j), j = 1, 2, ..., k. В этом состоит рекуррентный МНК. Для коррекции вектора BN не требуется ни обращения матриц, ни повторного формирования и решения нормальных уравнений. Метод может использоваться как при оперативной, так и при ретроспективной идентификации. Пусть 27 x N 1,1 x N 1,2 x .  x N 1, k Сформируем клеточные (или блочные) матрицы X ; xT X xT T  XT x ; Y y N 1 . Тогда  X  B N 1   XT x T  x   1 XT x y Y  ( X T X  xxT ) 1 ( XT Y  xy N 1 ) N 1 Воспользуемся далее формулой из леммы об обращении матриц (см. Современные методы идентификации систем. Под ред. П. Эйкхоффа. М.: Мир, 1983, стр. 391). Лемма. Если необходимые обратные матрицы существуют, тогда (A + BCD)1 = A1  A1B(C1 + DA1B) 1DA1 Доказательство. Если лемма справедлива, то по определению обратных матриц произведение исходной матрицы A + BCD на ее обратную должно дать единичную матрицу I: (A + BCD)[ A1  A1B(C1 + DA1B) 1DA1] = AA1 + BCDA1   AA1B(C1 + DA1B) 1DA1  BCDA1B(C1 + DA1B) 1DA1 = I + + BCDA1  B(C1 + DA1B) 1DA1  BCDA1B(C1 + DA1B) 1DA1 = = I+ BCDA1  BC[C1(C1 + DA1B) 1 + DA1B(C1 + DA1B) 1]DA1 = = I + BCDA1  BC[(C1 + DA1B)(C1 + DA1B) 1]DA1 = = I + BCDA1  BCDA1 = I, ч.т.д. Пусть теперь C = 1, B = b – матрица-столбец, D = bT, тогда можно записать ( A  bb T ) 1  A 1  A 1bb T A 1 1  b T A 1b . 28 Применяя последнее выражение, когда A = XTX; b = x, получим: B N 1  (PN   BN  PN xxT PN 1  x T PN x PN xxT B N 1  x T PN x )( XT Y  xy N 1 )   PN xy N 1   B N  PN x( y N 1  x T PN 1  x T PN x B N 1  B N  PN xxT PN 1  x T PN x xy N 1  xTB N 1  x T PN x PN x T xy N 1  1  x PN x ). Отсюда: ( y N 1  x T B N ). Т.о. для вычисления BN+1 необходимо знать BN; PN и новую информацию yN+1; x. Вычисления состоят только в умножении и сложении матриц, т.е. отсутствуют операции обращения или решения системы уравнений. Обратим внимание на то, что произведение x T B N  ŷ N 1, равно значению выходного сигнала, предсказываемое моделью (построенной за N шагов) по вектору новых значений x. Если это предсказываемое значение равно измеренному значению выхода, т.е. если xT BN = yN+1, то BN+1 = BN и никакого пересчета вектора BN не требуется. Для применения формулы пересчета вектора BN+1 на следующем (N+2)-ом шаге потребуется, очевидно, матрица PN+1. Эту матрицу можно получить коррекцией матрицы PN. Так как BN = PN XTY, то выражение можно записать как где BN+1 = (XTX + xxT)-1 (XTY + xyN+1) BN+1 = PN+1 (XTY + xyN+1), P xxT PN PN 1  ( X T X  xxT ) 1  PN  N . 1  x T PN x Последнее выражение, полученное с помощью леммы об обращении матриц, позволяет получить откорректированную матрицу PN. 29 19. Взвешенный МНК и другие разновидности МНК Взвешенный МНК применяется, когда различным ошибкам i придается различный вес или за счет старения результатов измерения от t1 к tN или за счет неодинаковой точности измерений. Метод состоит в использовании критерия вида: N J1   q i i2  E T QE; i 1 q1 где Q  0  q2      q i  0, 0  диагональная матрица весов.  qN J1  (Y  XB) T Q(Y  XB); B  (X T QX) 1 X T QY. Возможные принципы назначения весов: 1) Первым измерениям (наблюдениям) назначаются меньшие веса, последним большие (учет старения информации) при условии 0  qi  1. 2) Веса назначаются как величины пропорциональные обратной величине погрешности измерений: k k q i  ; или q i  , i i где i  инструментальная погрешность; i  среднеквадратичная погрешность измерений; k  коэффициент пропорциональности. Существуют и другие разновидности МНК Если в модели со многими входами ŷ( t i )  b1x1( t i )  b 2 x 2 ( t i )    b k x k ( t i ) (*) вместо xj(ti) взять сигнал одного (например, первого) входа в разных степенях, то можно построить методом МНК одномерную нелинейную модель ŷ( t i )  b1x1 ( t i )  b 2 x1 ( t i ) 2  b 3 x1 ( t i ) 3   Подставляя в модель (*) вместо xj(ti) различные комбинации x1(ti)p x2(ti)q ... xj(ti)r ..., можно строить МНК многомерные нелинейные модели. 30 В модели (*) вместо xj(ti) можно использовать различные функции от входных отсчетов xj(ti). Очень часто в этом случае используются системы ортогональных функций. 20. Получение модели по частотным характеристикам Моделирование с помощью активного эксперимента выполняется обычно вне контура регулирования. Наиболее распространены активные воздействия в виде синусоидальных сигналов, а также в виде ступенчатых и импульсных воздействий. Пусть на частотах i, i=1..N выполняется снятие частотных характеристик. xi = sin it W(s) yi = A(i) sin[it + (i)] Рис. 19. Схема снятия частотных характеристик По 2N результатам эксперимента A(i) = Ai; (i) = i вычислим значения вещественной и мнимой частотных характеристик Ri =Ai cos i; Ii = Ai sin i. Будем искать передаточную функцию в виде A 0  A1s    A ms m ; m  n. 1  B1s    Bn s n Соответствующее описание системы в виде дифференциального W (s)  dmx dn y dx dy    Bn n  A 0 x  A1    Am m . dt dt dt dt Пусть для краткости положим n = 3 и m = 2, тогда для каждого i-го уравнения: y  B1 опыта y = Ai sin(it + i); dy/dt = i Ai cos(it + i); d2y/dt2 = i2 Ai sin(it + i); d3y/dt3 = i3 Ai cos(it + i). x = sin it; dx/dt = i cos it; d2x/dt2 = i2 sin it; 31 После подстановки этих выражений для установившегося синусоидального режима в дифференциальное уравнение получим тождество, справедливое на данной частоте i для любого момента времени t: Ai sin(it + i) + B1 i Ai cos(it + i)  B2 i2 Ai sin(it + i)   B3 i3 Ai cos(it + i) = A0 sin it + A1 i cos it  A2 i2 sin it Учтем преобразования: Ai sin(it + i) = Ai (sin it cos i + cos it sin i) = Ri sin it + Ii cos it Ai cos(it + i) = Ai (cos it cos i  sin it sin i) = Ri cos it  Ii sin it Поскольку тождество имеет место для любого момента времени t, то выберем некоторый момент t* из условия: it* = /4, т.е. t* = /(4i), тогда sin it* = cos it* и тождество можно записать как Ri + Ii + B1 i (Ri  Ii)  B2 i2 (Ri + Ii)  B3 i3 (Ri  Ii) = A0 + A1 i  A2, или, собирая неизвестные в левой части: A0 + A1 i  A2  B1 i (Ri  Ii) + B2 i2 (Ri + Ii) + B3 i3 (Ri  Ii) = Ri + Ii. Если в полученную систему уравнений подставить экспериментальные значения Ri и Ii, то равенство для каждого i будет приближенным. Тогда систему уравнений можно решать МНК, введя невязку как разность между левой и правой частями. Введенная таким образом невязка будет линейна относительно искомых коэффициентов и вычисление их МНК даст несмещенные оценки. Точность вполне удовлетворительная (различие в коэффициентах обусловлено лишь погрешностями округления). 21. Идентификация систем по переходной характеристике Под переходной характеристикой h(t) понимается реакция динамической системы на единичный скачок. Преимуществом идентификации систем по переходной функции является то, что переходная функция легко может быть получена при любом 32 скачкообразном изменении входной величины на входе в пределах области линейности (переход с одного режима на другой). Часто переходную функцию можно получить из разгонной характеристики или из характеристики "выбега", т.е. из характеристик, получаемых при включении или выключении систем. Для достаточно быстродействующих систем переходную характеристику можно наблюдать на осциллографе при подаче на вход последовательности прямоугольных импульсов. Реальный скачок всегда имеет конечное время нарастания. Рис. 20. Функция реального единичного скачка Необходимо только стремиться к тому, чтобы время нарастания tнараст было бы значительно меньше, чем самая малая из постоянных времени, характеризующих систему: tнараст << Tmin. Т.е. источник скачка должен быть более быстродействующим. Идеальная единичная функция (идеальный скачок) имеет изображение по Лапласу: 1(t)  1/s. Соответствующая амплитудно-частотная характеристика Aид() = 1/  гипербола. Т.о. испытание динамической системы при помощи единичного скачка эквивалентно частотным испытаниям, когда основная энергия сигнала сосредоточена на низких частотах, а с повышением частоты энергия испытательного сигнала резко падает. 33 Рис. 21. АЧХ реального скачка. Пусть генератор скачка описывается как звено первого порядка с постоянной времени . Тогда  t ; 1 ; s(s  1) амплитудно-частотная характеристика выходного сигнала генератора: 1 Aреал()  X( j)  .  1  () 2 x (t )  1  e X(s)  При 0 Aреал()1/, т.е. Aреал() стремится к характеристике идеального скачка. При  Aреал()1/[()]. Таким образом, чтобы реальная характеристика приближалась к идеальной при повышении частоты необходимо уменьшать . Иначе говоря, чем более широкополосная система, тем более быстродействующим должен быть генератор испытательного сигнала. 22. Идентификация звена 1-го порядка по переходной функции Пусть система описывается дифференциальным уравнением первого порядка T y'(t) + y(t) = K x(t); W (s)  K ; Ts  1 34 при x(t) = 1(t) y( t )  K (1  e  t T )  h ( t ); lim h ( t )  K. t  1 h (T )  K (1  )  0.63K. e dh ( t ) K , т.е. касательная в начальной точке Производная t 0  dt T пересекает уровень К при t=T. Теоретически для определения параметра Т можно провести касательную в начальной точке и найти Т по точке пересечения касательной с уровнем К. Т.к. начальная часть h(t) обычно сильно зашумлена, то практически параметр Т с большей точностью определяется на уровне 0.63К. При t=T Рис. 22. Идентификация звена первого порядка. 23. Идентификация звена 1-го порядка с запаздыванием по переходной функции При наличии чистого запаздывания W (s)  K  s e ; Ts  1  t  h ( t )  K (1  e T ); t  ; вначале определяется время запаздывания по графику, а затем на уровне 0.63К определяется постоянная времени Т. 35 Рис. 23. Звено первого порядка с запаздыванием. Для определения Т и  можно использовать также уровни 0.33К и 0.7К (точнее: (1  e0.4)K и (1  e1.2)K). Если на уровне 0.33К имеем время t1, а на уровне 0.7 время t2, то параметры звена можно вычислить по формулам:  = 0.5 (3t1  t2); T = 1.25 (t2  t1). Рис. 24. Идентификация звена первого порядка с запаздыванием. Обозначим первый уровень a1K. а второй a2K, тогда   t1  t 2 a1K  K (1  e T ); a 2 K  K (1  e T ). Поставляя в эти выражения формулы (*), получаем: a1  1  e 0.5(3t1  t 2 )  t1 1.25( t 2  t1 )  1  e  0.4  0.3297; (*) 36 a2  1 e 0.5(3t1  t 2 )  t 2 1.25( t 2  t1 )  1  e 1.2  0.6988. 24. Идентификация параметров колебательного звена 2-го порядка Дифференциальное уравнение звена a2 y"(t) + a1 y'(t) + a0 y(t) = b0 x(t) . Для нулевых начальных условий: y(0) = 0; y'(0) = 0; звено можно описать передаточной функцией W ( p)  b0 Y ( p)  . X(p) a 0  a1p  a 2 p 2 Чаще эту передаточную функцию записывают в форме: W ( p)  K T02 p 2  2dT0 p  1  K02 p 2  2d0 p  02 , где K = b0/a0  статический коэффициент передачи; T0  1 / 0  a 2 / a 0  постоянная времени звена; 0  1 / T0  a 0 / a 2  собственная (резонансная) частота звена; d  a1 /(2 a 0a 2 )  степень затухания. Если известны три параметра K, T0, d, то коэффициенты b0, a1, a2 при произвольном значении a0 вычисляются по формулам b0 = K a0; a1 = 2d T0 a0; a2 = T02 a0. Характеристическое уравнение T02 p2 + 2d T0 p + 1 = 0. Корни характеристического уравнения p1,2   1 d2 d j    j, T0 T0 37 где  = d/T0  коэффициент затухания;   0 1 d 2  круговая частота колебаний (T  2 /   2T0 / 1  d 2  период колебаний); будут комплексными, т.е. звено будет колебательным, при 0  d < 1 . Переходная функция  K 1     h ( t )  L1   K 1  e t (cos t  sin t ) 2 2 p    1  2dT0 p  T0 p  изображена на рис. 25. Рис. 25. Измеряемые отрезки на колебательной переходной функции. [См. ТАУ, часть 1, под ред. А. В. Нетушила, 1976, стр.92] Для идентификации звена (т.е. для вычисления его параметров) по переходной функции h(t) определяются K, A1, A2, T и вычисляется затем степень затухания d в качестве идентифицируемых параметров по формуле: 1 d . 2   2 1    ln(A1 / A 2 )  одного из Действительно, рассмотрим выражение dh/dt = 0, из которого получим моменты времени, соответствующие экстремумам h(t). 38 d  {K[1  e t (cos t  sin t )]}  0. dt  Дифференцирование и сокращение дает (2 + 2) sin t = 0. Для первых двух максимумов h(t) имеем в результате t1m  ; t 2m  3. Подставляем t1m и t2m в h(t) и, вычитая из каждого максимума K, получим значения A1 и A2: A1  h ( t1m )  K  e   ; A 2  h ( t 2m )  K  e  3  Откуда ln(A1/A2) = 2/. Так как  /   d / 1  d 2 , то 2 d2  ln(A1 / A 2 )  ,    2   1  d2 A1 2d ;  A2 1  d2 1 . откуда d  2   2 1   ln( A / A ) 1 2   ln Кроме того, можно вычислить (2)  ln(A1 / A 2 ) 2 1 1  ;   ln(A1 / A 2 );   . T0 T T T Этими формулами можно пользоваться при малой степени затухания (примерно при d < 0.4). При больших d трудно определить A2, поэтому в этом случае определяют (см. рис. 25) K, A1, t1. 2 2 0  Тогда параметры могут быть вычислены как d 2 1 2 ; 0  ln(K / A1 ) 1  ln(K / A1 )  [  arctg ] 1   . t1        1   ln( K / A ) 1   Эти формулы применимы при 0.4 < d < 0.7. 39 25. Определение параметров апериодического звена 2-го порядка При коэффициенте затухания больше единицы звено второго порядка становится апериодическим и его передаточную функцию можно представить в виде: K W (s)  . (T1s  1)(T2 s  1) Соответствующая переходная функция: t t   T2 T1 e T2 ). e T1  h(t )  L    K (1  T2  T1 T1  T2  s  1  W (s)  Горизонтальная асимптота дает величину К. Т1 и Т2 находятся по точке перегиба и положению касательной в точке перегиба. Существует несколько методов определения T1 и Т2. Рис. 26. Апериодическое звено 2-го порядка Метод отрезков Tb и Tc Усовершенствованный метод идентификации апериодических объектов 2-го порядка (метод отрезков Tb и Tc). На экспериментальной h(t) (см. рис. 26) измеряется Тb и Тc = Т1 + Т2 и вводятся обозначения T TT T x  b ; y  1 22 ; d  1 . Tc T2 Tc Можно показать, что 40 d 1 d 1 d x d ; y . 1 d (1  d ) 2 По этим уравнениям можно построить кривую y = f(x) (для 1  x  e/2 и 0  y  0.25) и использовать ее для нахождения у по экспериментальному значению x. Рис. 27. Метод отрезков Tb и Tc.    Tc  T1  T2  y Решая систему уравнений T1T2 Tc2 относительно неизвестных T1 и T2, получаем: T1  Tc 1  1  4y 2  Tc d ; T2  Tc  T1. 1 d Если x > e/2  1.359, то система, возможно, не второго порядка. 26. Идентификация по апериодической переходной функции с точкой перегиба звена первого порядка с запаздыванием В некоторых случаях (например, при ориентировочных расчетах) можно h(t) апериодического звена второго или большего порядка аппроксимировать передаточной функцией вида: W (s)  K  s e . Ts  1 41 Проводя касательную в точке перегиба, как показано на рис. 26, в первом (грубом) приближении можно принять  = Ta; T = Tb (см. кривую 1 на рис. 28). Рис. 28. Аппроксимация с чистым запаздыванием Более подходящие значения  и Т могут быть найдены, если потребовать, чтобы аппроксимирующая ha(t) проходила через точку перегиба и чтобы касательная для h(t) в точке перегиба была бы также касательной и для ha(t). [См. ТАУ, часть 1, под ред. А.В.Нетушила, 1976, стр.263-264] При t > : h a ( t )  K (1  e  t T ). Положим h a (t 0 )  h 0 ; dh a ( t ) K , tt0  dt Tb тогда будем иметь систему уравнений:  t 0 K (1  e T K e T  t 0 T  )  h0    K    Tb  Решение системы: T  Tb (1  h0 ); K   t 0  Tln(1  h0 ). K 42 27. Метод кратных корней Определение постоянной времени Т и кратности n для модели системы в виде K W (s)  (Ts  1) n можно выполнить с использованием информации о всей кривой переходного процесса h(t). Метод использует вычисление площади, заключенной между уровнем установившегося значения К и кривой h(t):  S   [K  h ( t )]dt. Рис. 29. Площадь S. Поскольку теоретически для искомой модели переходная функция h(t) описывается выражением: t   n 1 t i   h ( t )  K 1  e T  , то i   i 0 T i!    t   t n 1 i t i t i n 11    t   t  n 11    t   i dt  K   e T   dt  KT   e T   d  T T T i  0 i! 0 i  0 i! 0 i  0 T i! В последнем выражении находим табличный интеграл: S  Ke T  x i  e x dx  (i  1)  i!. После подстановки табличного интеграла получаем: 43 S=KTn ~ Алгоритм вычисления T и n по экспериментальной кривой h ( t ) может быть следующим: ~ 1. По экспериментальной кривой h ( t ) или по статической характеристике оценивается величина К. 2. Каким-либо подходящим методом численного интегрирования вычисляется площадь ~ ~ t max S   [K  h ( t )[dt. 3. Для последовательных значений n = 1, 2, 3 ... вычисляют соответствующие значения T n 1,T n 2 ,T n 3 ,  по формуле ~ S T . nK 4. Для каждой пары значений n и T можно вычислить затем h(t) и ~ оценить тем или иным способом разность h(t)  h ( t ) , например:  1 t max t max ~ 2  [h ( t )  h ( t )] dt . 5. В качестве результата выбирается та пара n и T, которой соответствует минимум . 28. Метод площадей Метод ориентирован на вычисление характеристического уравнения для модели W (s)  K 1  A1s  A 2s 2  A 3s 3  коэффициентов , (*) Получение формул для вычисления коэффициентов Ai начнем с рассмотрения преобразования Лапласа: 1 A1  A 2s    h(t )  1 L 1   .   2 K s   s(1  A1s  A 2s  ) 1  A1s  A 2s 2   С другой стороны: 44  h ( t )    h ( t )  st L 1     1  e dt. K  0 K   Обозначим 1 - h(t)/K = x, тогда, учитывая, что e st  1   st (st ) 2 (st )3    , 1! 2! 3!   Lx   xdt  s  x  t ( t ) 2 dt  s 2  x dt    1! 2!  M 0  M1s  M 2s 2      M is i , i 0  ( t ) i dt. где M i   x i! Приравнивая результаты, получаем тождество: A1  A 2s  A3s 2   2 3  M 0  M1s  M 2s 2   1  A1s  A 2s  A3s   Приводя к общему знаменателю и приравнивая коэффициенты при равных степенях s, получаем: A1 = M0 A2 = M1 + A1M0 A3 = M2 + A1M1 + A2M0 A4 = M3 + A1M2 + A2M1 + A3M0                  An = Mn1 + A1Mn2 + … + An2M1 + An1M0; n = 2, 3, … Вычисляя по экспериментальной h(t) приближенные значения Mi, по полученным формулам можно последовательно вычислить Ai. Первые n уравнений в приведенной системе всегда имеют n неизвестных коэффициентов Ai. Если подставить в последнюю формулу выражение для Mi, то можно получить формулу для вычисления Ai в виде:  An   x     ( t ) n 1 ( t ) n  2 t dt  A1  x dt    A n  2  x dt  A n 1  xdt  (n  1)! 0 ( n  2)! 0 1!   x (( (( t t t t  A1 )  A2 )    A n 2 )  A n 1 )dt n 1 n2 n 3 1! Эта формула позволяет вычислять Ai по выражениям (см. лаб. раб.): 45  A1   xdt  A2   x(  t  A1 )dt 1 A3   x ((  t t  A1)  A 2 )dt 2 1 A 4   x ((( t t t  A1)  A2 )  A3 )dt 3 2 1 . . . . . . . . . . . . . . . . . . . . . . . .  t t t t A n   x (( ((  A1 )  A2 )    An 2 )  A n 1 )dt n 1 n2 n 3 1 29. Основное уравнение идентификации Часто подача активных воздействий на объект нежелательна. Тогда используют статистическую обработку сигналов, наблюдаемых на входе и выходе, и затем определяют характеристики объектов. Если входной сигнал x(t) = 0 при t < 0, то реакцию y(t) системы вычисляют с помощью интеграла свертки (свертка x(t) с весовой функцией g(t)): t t y( t )   g( t  ) x ()d   g() x ( t  )d. В общем случае, когда x(t) существует и при t < 0, необходимо вычислять y(t) по формуле y( t )  t  g( t  ) x ()d.  Выполним замену переменной интегрирования: положим  = t  z, тогда z = t  , d = dz, при  =  z = , при  = t z = 0 и   y( t )    g(z) x ( t  z)dz   g() x ( t  )d. Для стационарного эргодического случайного процесса с нулевым мат. ожиданием автокорреляционная функция 1 T  x ( t ) x ( t  )dt; T  2T T R xx ()  lim R xx ()  R xx (). (*) 46 Взаимная корреляционная функция двух процессов 1 T R yx ()  lim  y( t ) x ( t  )dt. T   2T  T Подставим сюда y(t)  1 T    g( ) x ( t  )d  x ( t  )dt.  T  2T  T  0  R yx ()  lim Изменим порядок интегрирования     1 T R yx   g()  lim x ( t  ) x ( t  )dt  d   g ()R xx (  )d.  T  2T  T  Последнее выражение записано с учетом того, что в формуле (*) аргумент Rxx есть разность аргументов функций под знаком интеграла. Учитывая второе из равенств (*), можно записать  R yx   g()R xx (  )d. или  R yx   g()R xx (  )d. Это соотношение носит название основное уравнение идентификации, или уравнение Винера-Хопфа. Пусть x(t) белый шум с интенсивностью с, тогда Rxx() = c(): R xx (  )  c(  );  R yx ()   g()c(  )d  cg() с учетом фильтрующего свойства -функции, т.е. с учетом того, что значение интеграла равно значению подинтегральной функции в точке, определяемой нулевым значением аргумента -функции. Таким образом, идентифицировать систему (т.е. вычислить, например, импульсную характеристику g(t)) можно просто наблюдая белый шум на входе и вычисляя взаимную корреляционную функцию выходного и входного сигналов. 30. Решение основного уравнения идентификации Уравнение Винера-Хопфа  R yx ()   g ()R xx (  )d относится к классу интегральных уравнений Фредгольма 1-го рода. 47 Конечно-разностный метод Верхний предел заменяется конечным T, т.е. решается уравнение T R yx ()   g ()R xx (  )d. Вводим дискретные значения i=i T, i=1...n; j=j T, j=1...n и записываем: n R yx ( j )   R xx ( j  i )g (i )T; j  1 n. i 1 В матричной форме Ryx = Rxx G, где R yx R yx (1 ) 1 R yx ( 2 ) ;  T  R yx ( n ) g(1 ) g G  (2 ) ;  g(n ) R xx (0) R xx (T)  R ( T ) R xx (0) xx R xx    R xx [(n  1)T ] R xx [(n  2)T ] Решение уравнения G = Rxx-1 Ryx.  R xx [(n  1)T ]  R xx [(n  2)T ] .   R xx (0)  Основная проблема: обращение матриц большой размерности (несколько сотен и более), что представляет трудности даже для обращения симметричных матриц. Кроме того, необходимо учесть, что Ryx и Rxx получаются из опыта по наблюдениям входа и выхода и потому дискретный аналог уравнения Винера-Хопфа следует записать в виде m y( j)   x ( j  i)g(i)   j ; j  1 n;x ( j  i)  x (i  j); i 1 или в матричной форме Y = X G + E, 1 где E   2 ; n  m.  n Тогда G можно найти из условия J = ET E = (Y - XG)T (Y - XG) = min. Положив как и ранее dJ/dG = 0 найдем G = (XTX)-1 XTY. 48 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. Вопросы к экзамену Два аспекта понятия моделирования. Понятие об идентификации. Причины необходимости создания новых моделей. Характеристики объектов, которые надо учитывать при создании моделей. Приемы упрощения моделей. Этапы построения моделей. Определение цели получения моделей. Определение ограничений и условий, учитываемых при построении моделей. Выбор подхода к решению задачи получения модели. Определение класса модели. Выбор метода решения задачи и ее решение. Общие принципы построения аналитических моделей. Модель поплавкового уровнемера. Модель процесса теплопередачи. Модель смесителя. Модель химического реактора. Модель емкости с изменяющимся уровнем. Метод наименьших квадратов в одномерном случае Метод наименьших квадратов в многомерном случае. Рекуррентный метод наименьших квадратов. Взвешенный метод наименьших квадратов и др. разновидности МНК. Получение модели по частотным характеристикам Свойства метода идентификации по частотным характеристикам (влияние различных факторов на погрешность идентификации). Идентификация систем по переходной характеристике Идентификация звена 1-го порядка по переходной функции Идентификация звена 1-го порядка с запаздыванием по переходной функции Идентификация параметров колебательного звена 2-го порядка Определение параметров апериодического звена 2-го порядка Идентификация звена 2-го порядка по переходной функции с запаздыванием Метод кратных корней Свойства кратных корней (влияние различных факторов на погрешность идентификации). Метод площадей Свойства метода площадей (влияние различных факторов на погрешность идентификации). Основное уравнение идентификации. Решение основного уравнения идентификации.