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

Клюенков Андрей Леонидович

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

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

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36565f1be77c40d58cfc

UUID: 31a56cbf-344c-4c10-b1c2-3b24c1928a07

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 10

Клюенков Андрей Леонидович

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 597812 bytes

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ КАФЕДРА МОДЕЛИРОВАНИЯ ЭКОНОМИЧЕСКИХ СИСТЕМ Клюенков Андрей Леонидович Выпускная квалификационная работа бакалавра Реализация алгоритмов оптимального управления с учетом реальных условий функционирования объектов управления Направление 010400 Прикладная математика и информатика Научный руководитель, доктор физ.-мат. наук, профессор Смирнов Н. В. Санкт-Петербург 2016

Содержание Введение. Обзор литературы 3 Постановка задачи 4 Глава 1. Решение задач оптимального управления при помощи адаптивного метода 6 §1.1. Сведение линейной непрерывной задачи оптимального управления к задаче линейного программирования . . . . . . . . . . . 6 §1.2. Сведение дискретной задачи оптимального управления к задаче линейного программирования . . . . . . . . . . . . . . . . . . . . 8 §1.3. Случай постоянных возмущений . . . . . . . . . . . . . . . . . . 9 §1.4. Идея адаптивного метода и его основные положения . . . . . . 10 §1.5. Адаптивный метод . . . . . . . . . . . . . . . . . . . . . . . . . . 12 §1.6. Построение начального плана . . . . . . . . . . . . . . . . . . . . 13 §1.7. Численная реализация . . . . . . . . . . . . . . . . . . . . . . . . 14 Глава 2. Прикладные задачи 16 §2.1. Стабилизация трехмассовой колебательной системы . . . . . . . 16 §2.2. Стабилизация системы при воздействии на нее постоянных возмущений . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 §2.3. Оптимальное управление макроэкономическими тенденциями на основе разностной схемы МОБ . . . . . . . . . . . . . . . . . 19 §2.4. Оптимальная стабилизация гибкого робота-манипулятора . . . 26 Заключение 31 Список литературы 32 Приложение 34 2

Введение. Обзор литературы В связи с возникновением необходимости решения различных технических и экономических задач, в середине прошлого века была создана теория оптимального управления. Направленная на оптимизацию функционалов, характеризующих всевозможные параметры математических моделей, данная теория позволяет при помощи теоретических подходов, а наряду с повсеместным внедрением в жизнь программных средств, и численно, находить наиболее выгодные режимы управления объектом. Под оптимальностью, как правило, понимается максимизация или минимизация некоторых характеристик объекта, таких как быстродействие, расход энергии, производительность. Для задач оптимального управления применяются различные динамические модели, описывающие поведение объектов управления. Примерами таких моделей служат системы обыкновенных дифференциальных уравнений, системы разностных уравнений, системы дифференциальных уравнений в частных производных. В XX веке различные ученые занимались развитием классической теории оптимального управления. Л. С. Понтрягин для широкого класса задач сформулировал принцип максимума [1]. Известным достижением в области оптимального управления стало динамическое программирование, разработанное Р. Беллманом [2]. Большой вклад в развитие теории внесли В. И. Зубов [3, 4, 5], Р. Е. Калман [6] и другие. С появлением программных средств разработки стало возможным создание регуляторов, позволяющих находить оптимальное управление с учетом реальных условий функционирования объектов управления. Для решения данных задач были предложены различные подходы и методы. Один из таких методов был разработан Р. Габасовым и его учениками [7, 8, 9]. Настоящая работа посвящена реализации алгоритмов оптимального управления в режиме реального времени. 3

Постановка задачи Рассмотрим систему дифференциальных уравнений общего вида ẋ = f (t, x, u), (1) x(t∗ ) = x0 , (2) с начальными условиями t ∈ [t∗ , t∗ ] = T, t∗ < t∗ < +∞, где x ∈ Rn , u ∈ Rr — вектор управляющих параметров, f (t, x, u) — непрерывная по всем аргументам и непрерывно-дифференцируемая по x и u векторфункция. Множество всех допустимых управлений обозначим U . Функция u(t) ∈ U, t ∈ T , называется дискретным управлением (с периодом квантования h), если  u(t) = u(t∗ + kh) = ukh , t ∈ [t∗ + kh, t∗ + (k + 1)h[ k = 0, N − 1 (3) где h = (t∗ − t∗ )/N — период квантования, N — натуральное число. И кроме того, l∗ 6 u(t) 6 l∗ , при t ∈ T . Множество точек разбиения будем обозначать Tu = {t∗ , t∗ + h, . . . , t∗ − h}. Сформулируем задачу терминального управления. Для решения системы (1) с начальными условиями (2), в классе кусочно-постоянных управлений (3), обеспечить максимум целевой функции cT x(t∗ ) → max, (4) Hx(t∗ ) = g, (5) при условии где g ∈ Rm , rank H = m 6 n — в общем случае. Таким образом уравнение (5) задает линейное многообразие в пространстве Rm , на которое необходимо вывести программное движение системы (1) в конечный момент t = t∗ . 4

Помимо задач оптимального управления с функционалом (4) будем рассматривать задачи с функционалами непосредственно описывающими качество управления: Z t∗ cT u(t)dt → max . t∗ Далее в работе будет рассмотрен метод решения поставленной задачи оптимального управления и некоторые его приложения. 5

Глава 1. Решение задач оптимального управления при помощи адаптивного метода Для поиска оптимального управления проведем сведение задачи оптимального управления к задаче линейного программирования, которую, в свою очередь, будем решать при помощи адаптивного метода [7]. §1.1. Сведение линейной непрерывной задачи оптимального управления к задаче линейного программирования Рассмотрим систему ẋ = A(t)x + b(t)u, (6) с начальными условиями (2). Запишем общее решение системы (6), (2) в форме Коши в момент t = t∗ x(t∗ , t∗ , x0 ) = Y (t∗ )[Y −1 (t∗ )x0 + Zt∗ Y −1 (τ )b(τ )u(τ )dτ ], (7) t∗ где Y (t) — фундаментальная матрица соответствующей однородной системы. Тогда, подставляя полученное представление (7) в функционал (4), получим эквивалентное представление Zt∗ cT Y (t∗ )Y −1 (τ )b(τ )u(τ )dτ → max . (8) t∗ Рассмотрим сопряженную систему с начальными условиями ż = −AT (t)z, (9) z(t∗ ) = c, где c и t∗ из условий (4), (5). Ее решение представимо в виде zc (t) = Z(t)Z −1 (t∗ )c, где Z(t) — фундаментальная матрица системы (9). Запишем основное свойство этой фундаментальной матрицы: Z T (t) = Y −1 (t). 6

Используя данное свойство, легко убедиться, что zcT (t) = cT (Z −1 (t∗ ))T Z T (t) = cT Y (t∗ )Y −1 (t). (10) Учитывая (10), условие (8) можно записать в виде Zt∗ zcT (τ )b(τ )u(τ )dτ → max . (11) t∗ Таким образом, мы можем не вычислять фундаментальную матрицу системы (6), а найти решение сопряженной системы численным методом. Обозначая строки матрицы H за hi (i = 1, m) и проводя рассуждения аналогичные рассуждениям из предыдущего абзаца, можем найти численное решение для ограничения (5), которое с учетом (7) может быть записано в виде Zt∗ HY (t∗ )Y −1 (τ )b(τ )u(τ )dτ = g0 , (12) t∗ где g0 = g − HY (t∗ )Y −1 (t∗ )x0 — m-мерный вектор. Учитывая, что мы ищем управление в виде (3), можем записать (11) и (12) соответственно следующим образом X γk uk → max, (13) k∈Tu X dk uk = g0 , (14) k∈Tu где dk = t+h R HY (t∗ )Y −1 (τ )b(τ )dτ , γk = t t+h R zcT (τ )b(τ )dτ , при t = t∗ + kh, t k = 0, N − 1. Запишем полученные выражения в векторной форме. Для этого введем обозначения: ω = (u0h , u1h , . . . , ukh , . . . , u(N −1)h )T , C = (γ0 , γ1 , . . . , γN −1 ), D = (d0 , d1 , . . . , dN −1 ). Тогда окончательно задачу (13)–(14) можно записать 7

в виде    Cω → max,    Dω = g0 ,     L∗ 6 uk 6 L∗ , k = 1, N , (15) где L∗ = (l∗ , . . . , l∗ ), а L∗ = (l∗ , . . . , l∗ ). Таким образом, мы свели задачу оптимального управления к желаемому виду (15) — канонической задаче оптимального управления. §1.2. Сведение дискретной задачи оптимального управления к задаче линейного программирования Рассмотрим разностную управляемую систему x(t + 1) = Ax(t) + bu(t), (16) x(0) = x0 , (17) с начальными условиями здесь x, b ∈ Rn , u ∈ R1 , A — (n×n)-матрица. Будем предполагать, что время t может принимать следующие значения t ∈ {0, 1, . . . , N } = T0 . Поставим задачу терминального управления. Обеспечить максимум функционала cT x(N ) → max (18) Hx(N ) = g, (19) с ограничениями где rank H = m 6 n, g ∈ Rm . Управление будем считать ограниченным |u(t)| 6 L для всех t ∈ T0 . Введем также в рассмотрение множество T1 = {0, 1, . . . , N − 1}. Выпишем общее решение задачи (16)–(17) t x(t) = A x0 + t X i=1 8 Ai−1 bu(t − i). (20)

В конечный момент времени t = N решение (20) будет иметь вид N x(N ) = A x0 + N X Ai−1 bu(N − i). (21) i=1 Подставляя (21) в (18) N −1 X cT Ai bu(N − i) → max . i=0 Теперь подставим (21) в ограничения (19), получим N −1 X HAi bu(N − i) = g0 , i=0 где g0 = g − HAN x0 . Введем обозначения γt = cT At b, dt = HAt b, тогда X γt u(N − t) → max, t∈T1 X dt u(N − t) = g0 . t∈T1 Тогда полагая ω = (u(N − 1), u(N − 2), . . . , u(0))T , C = (γ0 , γ1 , . . . , γN −1 ), D = (d0 , d1 , . . . , dN −1 ), придем к виду (15). §1.3. Случай постоянных возмущений Рассмотрим теперь систему линейных дифференциальных уравнений в условиях постоянно действующих периодических возмущений f (t) ẋ = A(t)x + b(t)u + f (t), (22) x(τ ) = z, Hx(t∗ ) = g, 0 6 u(t) 6 L, t ∈ T (τ ) = [τ, t∗ ], (23) где x ∈ Rn , A(t) — (n × n)-матрица, b(t) — n-мерная вектор-функция, u ∈ R1 — управляющий параметр. Будем искать управление в классе кусочнопостоянных функций, а множество допустимых управлений обозначим U (τ ). Пусть u0 (t|(τ, z)) ∈ U (τ ), t ∈ T (τ ) — оптимальное управление задачи (22), 9

(23) для (t, z), z ∈ X(τ ), X(τ ) — множество состояний, для которых (22), (23) имеет решение при фиксированном τ . Предположим, что в моменты времени t∗ , t∗ +h, . . . , τ определены управления u∗ (t∗ ), u∗ (t∗ + h), . . . , u∗ (τ ). Тогда реальное состояние системы x(τ + h), в которое она перешла за время h, связано с состоянием в отсутствии возмущения x0 (τ + h) в момент τ + h следующим образом: Z τ +h ∗ 0 x (τ + h) = x (τ + h) + Y (τ + h)Y −1 (s)f (s) ds. τ В момент τ + h надо определить оптимальное управление: u∗ (τ + h) = u0 (τ + h|(τ + h, x∗ (τ + h)), cT x(t∗ ) → max, ẋ = A(t)x + b(t)u, x(τ + h) = x∗ (τ + h), Hx(t∗ ) = g, t ∈ T (τ + h). 0 6 u(t) 6 L, Таким образом, в момент времени τ +h надо решить следующую задачу линейного программирования N −1 X ck u(k) → max, k=τ N −1 X dk u(k) = g0 (τ ) + ∆g(τ ), k=τ ∆g(τ ) = Y (t∗ )Y −1 (τ )x∗ (τ ) + dk (τ )u∗ (τ ) − Y (t∗ )Y −1 (τ + h)x∗ (τ + h). Полученная задача линейного программирования может быть решена при помощи адаптивного метода. §1.4. Идея адаптивного метода и его основные положения Обозначим через I множество индексов строк матрицы D, а через J — множество индексов ее столбцов. Выделим подмножества Ibs ∈ I и Jbs ∈ J с одинаковым количеством элементов в каждом из них (|Ibs | = |Jbs |) и, на основе этих подмножеств, составим квадратную матрицу Dbs = D(Ibs , Jbs ). Совокупность множеств Kbs = {Ibs , Jbs } будем называть опорой, а матрицу 10

Dbs — опорной (базисной) матрицей, если det Dbs 6= 0. Также введем неопорные элементы: In = I \ Ibs , Jn =! J \ Jbs , а также представим матрицу D Dbs Db,n в блочном виде D = . Dn,b Dn Будем рассматривать задачу линейного программирования (15). Пусть ω — базисный план этой задачи, Dbs — его базисная матрица, ω — другой базисный план, ∆ω = ω − ω — приращение плана ω. Подсчитаем приращение целевой функции ∆ϕ = Cω − Cω = C∆ω. (24) Учитывая, что ω и ω— планы задачи (15), D∆ω = Dω − Dω = b − b = 0 или в компонентной форме Dbs ∆ωbs + Dn ∆ωn = 0, откуда следует, что −1 ∆ωbs = Dbs Dn ωn . (25) Подставим теперь значение приращения по опорным компонентам вектора ω из уравнения (25) в равенство (24) −1 −1 C∆ω = Cbs ∆ωbs + Cn ∆ωn = Cbs Dbs ∆zbs + (Cn − Cbs Dbs Db,n )∆ωn . Вводя обозначения −1 vbs = Cbs Dbs , ∆n = Cn − vbs Db,n , получим формулу приращения целевой функции X C∆ω = ∆n ∆ωn = ∆j ∆ωj . j∈Jn Физический смысл оценок и потенциалов подробно рассмотрен в работе [8]. Таким образом можно сформулировать следующую теорему [9]. Теорема 1. Для оптимальности плана ω необходимо и достаточно 11

существования такой опоры Kbs , при которой выполняется соотношение  ∆j 6 0 при ωj = l∗ , j ∈ Jn . ∆ > 0 при ω = l∗ , j j Адаптивный метод является итеративным. Последовательность приближений можно описать в виде следующей схемы 1 2 n {ω 1 , Kbs } → {ω 2 , Kbs } → · · · → {ω n , Kbs }. 1 } первой фазой, а Будем называть построение первого приближения {ω 1 , Kbs последовательности всех остальных — второй. §1.5. Адаптивный метод Приближения базисного плана можно представить в виде ω = ω + ∆ω, где ω — новый базисный план. Учитывая идею метода, приращение ∆ω можно представить в виде ∆ω = θ0 l. (26) Вектор l ∈ Rn называется направлением изменения плана ω, θ0 > θ — шаг вдоль этого направления. Из свойства D∆ω = 0 следует, что Dbs lbs + Dn ln = 0, а тогда для направления по базисным переменным справедлива формула −1 lbs = −Dbs Dn ln . (27) Из соотношений (26) и (27) видно, что для любого выбора ln основные ограничения задачи для ω выполняются, тогда не умаляя общности можем выбрать вектор неопорных компонент вектора l следующим образом ln = (0, . . . , 0, lj0 , 0, . . . , 0), где для компоненты lj0 , находящейся на позиции j0 ∈ Jn , выполняется lj0 = sign∆j0 . 12

При вычислении оптимальной длины шага θ0 разумно исходить из соображений о том, что вдоль направления l целевая функция возрастает. Таким образом, будем увеличивать значение θ до тех пор, пока не будет достигнута граница допустимых значений множества планов. Такие требования могут быть достигнуты при выборе θ0 = min{θj0 , θj∗ }, где θj∗ = min θj , θj0 = l∗ − l∗ . j∈Jbs Для завершения итерации требуется определить новое опорное множество, а также проверить условие оптимальности построенного плана. Перейдем к рассмотрению этих вопросов. Важным в вопросе выбора опоры является выбор длины шага θ0 . В случае, если оптимальная длина шага равна θj∗ , следует заменить индекс j∗ из опорного множества, при котором новый план достигает границы допустимого множества, на индекс j0 . Тогда новое множество опорных индексов будет иметь вид J bs = (Jbs \ j∗ ) ∪ j0 . Когда же при выборе длины шага θ0 активным становится ограничение j0 , то есть нарушаются ограничения при неопорных компонентах плана, новое множество неопорных компонент примет вид J n = Jn \j0 . Если такое множество оказывается пустым J n = ∅, это говорит о том, что новый план является оптимальным. §1.6. Построение начального плана Для того чтобы приступить к поиску оптимального плана задачи (15), необходимо определить начальный план задачи, который будет являться допустимым. Стоит отметить, что адаптивный метод позволяет использовать в качестве начального плана не только вычисленные на первой фазе значения, но и использовать апостериорную информацию, полученную от специалистов, так как данный метод оперирует не только вершинами симплекса, но и точками, находящимися внутри него. Целью первой фазы является поиск начального базисного плана или, в 13

случае, если такого плана не существует, доказательство несовместности ограничений задачи. Разобьем множество индексов J на два подмножества J∗ и J ∗ таким образом, чтобы J∗ ∩ J ∗ = ∅, J∗ ∪ J ∗ = J и положим ω̃ = (ω̃j = l∗ , j ∈ J ∗ ; ω̃j = l∗ , j ∈ J∗ ). Вычислим теперь вектор невязок ξ = g0 − Dω̃. Обозначим I + = {i ∈ I : ξ > 0}, I − = {i ∈ I : ξ < 0}. Если множество I − = ∅, то задача первой фазы решена, а вектор (ω̃, ω̃f = g0 − Dω̃) — базисный план задачи с базисной матрицей Dbs = (an+i , i ∈ I) = E. Предположим, что I − 6= ∅. Введем вспомогательную задачу для первой фазы X ωn+i → max, − i∈I −  D(i, J)ω + ωn+i = g0i , i ∈ I + , D(i, J)ω − ω − n+i = g0i , i ∈ I , (28) L∗ 6 ω 6 L∗ , 0 6 ωn+i 6 l∗ , i ∈ I + , 0 6 ωn+i 6 |ξi |, i ∈ I − , где ωn+i , i ∈ I + — свободные переменные, а ωn+i , i ∈ I + — искусственные переменные. Вектор (ω̃, ωn+i = ξi , i ∈ I + , ωn+i |ξi |, i ∈ I − ) является базисным планом задачи (28) с базисной матрицей D̃bs = (an+i = ei , i ∈ I + , an+i = = −ei , i ∈ I − ). Взяв его в качестве начального плана и применив к задаче (28) вторую фазу метода, можем построить оптимальный базисный план ∗ ∗ (ω ∗ , ωn+i , i ∈ I) с оптимальной базисной матрицей Dbs = (aj , j ∈ Jbs ). §1.7 Численная реализация Основываясь на данном алгоритме, был написан комплекс программ в пакете MATLAB. Для поиска численного значения определенного интеграла используется метод Симпсона, а для решения системы линейных дифферецниальных уравнений метод Рунге — Кутты четвертого порядка точности. 14

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

Глава 2. Прикладные задачи §2.1. Стабилизация трехмассовой колебательной системы В качестве примера, иллюстрирующего u x1 работу алгоритма для линейных систем, расm1 смотрим задачу стабилизации колебательc1 ной трехмассовой системы за фиксирован- x2 m2 ное время с минимальным расходом топлива c2 (см. рис. 1). x3 m3 Математическая модель задачи имеет вид c3 Z25 Рис. 1. Трехмассовая колебательная система u(t)dt → min, 0 ẋ1 = x4 , ẋ2 = x5 , ẋ3 = x6 ẋ4 = −2,33x1 − x2 + 0,33x3 + u, ẋ5 = −0,88x1 − 2,06x2 − 0,88x3 , x1 (0) = x2 (0) = x3 (0) = 0, ẋ6 = 0,25x1 − 0,75x2 − 1,75x3 , x4 (0) = 2, x5 (0) = 1, x6 (0) = 1, x1 (25) = x2 (25) = x3 (25) = x4 (25) = x5 (25) = x6 (25) = 0, 0 6 u(t) 6 1, t ∈ [0, 25[, где x1 = x1 (t), x2 = x2 (t), x3 = x3 (t) — отклонение от положения равнове1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 25 Рис. 2. Оптимальное управление сия первой, второй и третьей масс соответственно, u = u(t) — управляющий параметр силы. 16

Для построения оптимального управления был использован адаптивный метод. Результаты вычислений представлены на рис. 2–3. 2 1.5 1 0.5 0 −0.5 −1 −1.5 −2 0 5 10 15 20 25 Рис. 3. Интегральные кривые замкнутой системы §2.2. Стабилизация двухмассовой колебательной системы при воздействии на нее постоянных возмущений Рассмотрим задачу стабилизаu ции с минимальным расходом топливом колебательной двухмассовой си- x1 m1 стемы (рис.4), на которую оказывает влияние постоянно действующее воз- x2 мущение f (t) = 0,2 sin 10t. m2 Такую колебателную систему можно описать в виде соотношений Рис. 4. Двухмассовая колебательная система 17

Z 20 u(t)dt → min, ẋ1 = x3 , ẋ2 = x4 , 0 ẋ3 = −x1 + x2 + u, ẋ4 = 0,1x1 − 1,02x2 + 0,2 sin 10t, x(0) = (0; 0; 2; 1), x(20) = (0; 0; 0; 0), 0 ≤ u(t) ≤ 1, t ∈ [0, 20]. Данная задача была решена предложенным методом, с учетом действующего возмущения, и проведено сравнение с моделью без возмущений. Полученные результаты представлены на рисунках 5–7. 1 0.8 0.6 0.4 0.2 0 0 2 4 6 8 10 12 14 16 18 20 14 16 18 20 Рис. 5. Оптимальное управление 3 2.5 2 1.5 1 0.5 0 −0.5 −1 −1.5 −2 0 2 4 6 8 10 12 Рис. 6. Интегральные кривые замкнутой системы На рис. 7 показан фазовый портрет системы для переменной, на которую действует возмущение (здесь x2 — положение второй массы, а x4 — ее скорость). 18

1.5 1 0.5 0 −0.5 −1 −1 −0.5 0 0.5 1 1.5 Рис. 7. Фазовый портрет системы по переменным x2 , x4 (сплошной линией изображена кривая для возмущенной, пунктирной – без возмущений) §2.3. Оптимальное управление макроэкономическими тенденциями на основе разностной схемы МОБ Рассмотрим применение адаптивного метода для модели межотраслевого баланса (МОБ). Таблица МОБ и обозначения. Таблица МОБ описывает n секторов экономики региона. Рассмотрим таблицу межотраслевого баланса, которая описывает n секторов региональной экономики. Таблица межотраслевого баланса состоит из 4 квадрантов (см. таблицу 1) [10]. Рассмотрим каждый из них более подробно. Будем использовать обозначения в соответствии с работой [11]. Первый квадрант представляет собой (n×n)-матрицу производственной сферы Ap = {pij }. Столбцы определяют промежуточное потребление каждой j-ой отрасли как производителя. Элементы матрицы можно представить в следующем виде pij = Pi aij Inj , где Pi – это цены потребленной продукции каждой i -ой отрасли, aij – технологический коэффициент и, наконец, Inj – 19

объём годового выпуска j-ой отрасли в натуральном выражении. Производящие Потребляющие отрасли Конечный Валовой отрасли 1 2 3 ... n продукт (Y) продукт (X) 1 p11 p12 p13 ... p1n Y1 X1 2 p21 p22 p23 ... p2n Y2 X2 3 p31 p32 p33 ... p3n Y3 X3 . . . . . . . . . . . . (II) . . . . . . . . n pn1 pn2 pn3 ... pnn Yn Xn Оплата труда W1 W2 W3 ... Wn Оплата налогов T x1 T x2 T x3 (III) T xn Чистая прибыль P rh1 P rh2 P rh3 ... P rhn Валовой продукт X1 X2 X3 ... Xn (I) Vb (IV ) Таблица 1. Таблица межотраслевого баланса Рассмотрим суммы по столбцам матрицы Ap. Они равны стоимости промежуточного потребления в j-ой отрасли n X P pj = Inj Pi aij . i=1 Сумма по каждой i-ой строке равна стоимости объема реализованной продукции Xi этой отрасли экономики X i = Pi n X ! (29) aij Inj + Y nj j=1 Второй квадрант матрицы МОБ — (n + 1)-ый столбец, представляет собой n-мерный вектор-столбец стоимостей продукции конечного потребления Y = (Y1 , . . . , Yn )T = (P1 Y n1 , . . . , Pn Y nn )T . Третьим квадрантом матрицы «затраты-выпуск» является (n + 1)-ая строка V. Ее элементы Vj представляют собой показатели добавленной стоимости Vj = Ij − P pj = Inj Pj − n X ! aij Pi . i=1 Четвертым квадрантом матрицы МОБ является государственный бюджет Vb , который представлен в (n + 1)-ом диагональном элементе матрицы 20

МОБ. Он формируется как сумма всех налогов и других выплат. Как видно из приведенной таблицы, каждый элемент строки «Итого» (Валовой продукт) равен сумме затрат на промежуточное потребление и добавленную стоимость. Он представляет собой ожидаемую стоимость выпущенной продукции в j-ой отрасли, т. е. Ij = P pj + Vj . Тогда, учитывая соотношения для выпусков, получим Ii = Pi Ini = Pi n X ! aij Inj + Y ni . j=1 Таким образом, можно заметить, что V = Y . Основное уравнение баланса в натуральной форме (E − A)In = Yn. Для построения балансовых уравнений введем в рассмотрение матрицу относительных цен R{rij } rij = pij , pij = Pi aij Inj . Ij Тогда система уравнений (29), для случая (Ii 6= Xi ), в новых обозначениях примет следующий вид    r11 I1 + r12 I2 + · · · + r1n In + Y1 = X1 ,      r21 I1 + r22 I2 + · · · + r2n In + Y2 = X2 ,   ...      rn1 I1 + rn2 I2 + · · · + rnn In + Yn = Xn , или в матричной форме RI + Y = X. (30) В случае динамического равновесия (Ii = Xi ) система уравнений (30) принимает следующий вид (E − R)I = Y, 21

где I = (P1 In1 , . . . , Pn Inn ), Y = (P1 Y n1 , . . . , Pn Y nn ). Полученное уравнение является основным уравнением баланса в стоимостной форме. Разностная модель МОБ Предположим, что годовой объем продаж равен объему годового выпуска продукции X = I, а общая добавленная стоимость равна общему потреблению V = Y. Ii = ri1 I1 + ri2 I2 + · · · + rin In + Y ri In+1 . (31) В уравнении (31) rij = pij /Ij это элементы матрицы относительных цен R = {rij }, Y ri = Yi /In+1 итоговое потребление, нормированное по ВВП. Под ВВП мы понимаем сумму добавленной стоимости Vj и бюджета Vb , который является добавленной стоимостью сферы потребления. Чистая операционная прибыль P rhj – это источник инвестиций Cpj = P rhj . Значение P rhj отличается от добавленной стоимости Vj , оплаты труда работника Wj и налогов T xj . Обозначим изменение объема выпуска продукции за единицу времени как ∆Ij = Ij (t) − Ij (t − 1). Инвестиции необходимые для расширения производства (ускорения) выпуска продукции можно считать пропорциональными желаемому ускорению. Таким образом получим систему разностных уравнений уравнений P rhj = F ej ∆Ij , j = 1, . . . , n. (32) где F ej представляют капитальные фондоемкости каждого сектора. Они являются коэффициентами пропорциональности между ростом выпуска продукции и необходимыми для этого инвестициями. Уравнения, аналогичные (32), можно записать для In+1 (ВВП): P rhb = F eb ∆In+1 . (33) Здесь F eb – это фондоемкость сферы потребления. Для вычисления F ej можно использовать следующие выражения: P rhj (T2 )(T2 − T1 ) . F ej = Ij (T2 ) − Ij (T1 ) 22

Используя (31)–(33) и вводя обозначение I = (I1 , . . . , In , In+1 ), получаем систему уравнений e + E, I(t) = DI(t − 1), D = MR где E – единичная матрица,  r11 r12 · · ·   r21 r22   .. e R =  ... .   r  n1 rn2 V r1 V r1      M=     r1n Y r1 r2n Y r2 .. .. . . ··· rnn  V r1 − W r1 ··· F e1 .. ... . 0 0      ,  Y rn   rg ··· ... · · · V rn (34) 0 .. . V rn − W rn ··· F en ··· 0 0 .. .     .  0   V rb − W rb  F eb Здесь rg обозначает ставку обобщенного налога, которая определяется как доля бюджета ВВП, V rj = Vj /Ij и W rj = Wj /Ij это относительные добавленные стоимости и относительные оплаты труда работников. В систему (34) может быть введено управление как экзогенный параметр. В частности, таким параметром могут выступать инвестиции. Тогда динамику системы можно описать следующими уравнениями I(t) = DI(t − 1) + QU(t − 1), (35) 0 ≤ uj ≤ Lj . (36) Здесь U = (u1 , . . . , un+1 )T — вектор управления (инвестиции), Lj > 0 — константы, элементы матрицы Q = {qij }n+1 i,j=1 , где qij = δij /F ej и δij — это дельты Кронекера, i, j = 1, . . . , n + 1. 23

Прогнозирование и сравнение Стоит отметить, что для полученной модели МОБ, решения будут отличаться от решений для модели, описываемой системой дифференциальных уравнений [12]. В связи с этим проведем сравнение этих двух моделей. в млн руб. 1 2 3 Y 1 9 069 545 1 774 978 3 740 019 10 938 888 2 716 976 633 869 209 075 6 476 850 3 1 615 173 893 268 4 150 585 6 074 466 V 7 538 260 4 171 580 11 825 365 6 170 500 W 1 975 907 2 168 750 3 838 591 1 057 900 Таблица 2. Агрегированная таблица за 2006 год Рассмотрим таблицу межотраслевого баланса за 2006 год. Для этих целей воспользуемся агрегированными данными по трем секторам (промышленность, жизнеобеспечение и инфраструктура), представленные Россстатом (табл. 2) [13]. Для проверки адекватности модели ее матрица D построена по таблице товаров и услуг за 2005 год. Фондоемкости вычислены с использованием данных по выпуску и ВВП за 2004 и 2005 годы: F ej (2005) = P rhj (2005) . Ij (2005) − Ij (2004) Результаты прогнозирования представлены в таблице 3 I(2004) I(2005) I(2006) bI(2006) eI(2006) 11 994 582 15 547 550 18 939 954 20 790 352 19 984 787 5 126 942 6 067 326 7 473 695 7 218 155 7 599 113 12 415 833 15 478 184 19 925 044 17 704 991 18 393 854 18 611 055 24 227 335 29 705 705 30 452 286 29 720 487 Таблица 3. Реальные и смоделированные выпуски в текущих ценах Здесь bI(2006) – спрогнозированные значения выпусков для системы дифференциальных уравнений, а eI(2006) – для разностных. Для сравнения результатов прогнозирования с реальными данными за 24

2006 год использовалась формула относительной ошибки ||bI(2006) − I(2006)|| r= 100%. ||bI(2006)|| В таблице 4 используются следующие обозначения: r1 – относительная ошибка прогнозирования дифференциальной системы, а r2 – для разностной. Таким образом можно увидеть, что разностная модель наиболее точно описывает динамику выпусков продукции на данном промежутке времени. r1 r2 7,28 % 4,53 % Таблица 4. Относительная ошибка прогнозирования Оптимальное управление инвестициями в экономике Поставим задачу оптимального управления, которая заключается в увеличении выпусков и ВВП на 8% за год при минимальных инвестициях 4 X (u1 (t) + u2 (t) + u3 (t) + u4 (t)) → min, uj t=1 I(t) = D2006 I(t − 1) + Q2006 U(t − 1), 4 X 4 X uj (t) ≤ 5 · 1011 , I(0) = I2006 , I(1) = I2007 = 1.08I(0), t=1 j=1 где матрицы D2006 , Q2006 , I2006 получены из таблиц 1, 2. Предполагается, что инвестиции поступают в экономику ежеквартально, а их сумма не превосходит 0,5 триллионов рублей в квартал. Решение данной задачи найдено при помощи адаптивного метода, результаты представлены в таблице 5. u∗1 u∗2 u∗3 u∗4 1-ый квартал 2-ой квартал 3-ий квартал 4-ый квартал 0 0 0 2.45 ∗ 1011 0 0 0 1.61 ∗ 1011 0 0 4.61 ∗ 1011 0.87 ∗ 1011 5 ∗ 1011 5 ∗ 1011 0.39 ∗ 1011 0 Таблица 5. Оптимальное распределение инвестиций 25

Функционал оптимальности при полученном распределении инвестиций имеет значение 4 X 4 X uj (t) = 1.99 ∗ 1012 . t=1 j=1 §2.4. Оптимальная стабилизация гибкого робота-манипулятора Перейдем теперь к вопросу о построении оптимального управления для нелинейных систем дифференциальных уравнений. В общем случае для систем, в которых нелинейность по фазовым переменным выражается в виде периодических функций или функций, которые можно заменить их кусочнолинейной аппроксимацией, также можно применять предложенный алгоритм. Рассмотрим модель односвязного робота-манипулятора с парой моторов [14, c. 313]. Когда соединения рассматриваются как пружины кручения, система динамических уравнений принимает вид φ2 + mgd cos φ1 = 0, J1 φ̈1 + F1 φ̇1 + K φ1 − N φ2 di K J2 φ̈2 + F2 φ̇2 − φ1 − = Kt i, L + Ri + Kb φ̇2 = u, N N dt (37) где φ1 и φ2 — углы поворота соединений на каждом из моторов, i — ток якоря, u — напряжение якоря, R и L — сопротивление якоря и его индукция, J1 и J2 — моменты инерции, F1 , F2 — коэффициенты вязкого трения, K — коэффициент жесткости пружины, Kt — крутящий момент, Kb — коэффициент обратной ЭДС, m — масса соединения, d — положение центра тяжести соединения, N — передаточное число и g — ускорение свободного падения. Преобразуя систему (37), как это показано в [14], получим ẋ1 = x2 + θ1 x1 , ẋ2 = x3 + θ2 x1 + θ3 cos x1 , ẋ3 = x4 + θ4 x1 + θ5 cos x1 , ẋ4 = x5 + θ6 x1 + θ7 cos x1 , ẋ5 = θ8 cos x1 + b0 u, 26 (38)

где R F1 F 2 + + θ1 = − L J1 J2 Kb Kt K F1 F2 K R F1 F 2 + + + + + θ2 = − L J1 J2 L J2 J1 J2 N 2 J1 J2 mgd θ3 = − J1 K F1 F2 F1 K F2 K K b F 1 K 1 R K + + + + + θ4 = − L J1 J2 N 2 J1 J2 J1 J2 N 2 J1 J2 L J1 J2 mgd R F2 + θ5 = − J1 L J2 F1 K F2 K Kb F1 K1 R + + θ6 = − L J1 J2 N 2 J1 J2 L J1 J2 K RF2 Kb Kt mgd θ7 = − + + N2 L L J1 J2 R mgdK Kt K θ8 = − , b0 = > 0. 2 L J1 J2 N J1 J2 N L Для (38) рассмотрим задачу стабилизации с минимальным расходом топлива B(t∗ ) = Zt∗ u(t) dt → min 0 и краевыми значениями x(0) = z, x(t∗ ) = 0, где x = (x1 , . . . , x5 )T , z = (z1 , . . . , z5 )T , t ∈ T = [0, t∗ ]. Кроме того, будем полагать управление ограниченным |u(t)| < L = const. Заменим нелинейную характеристику системы cos x на периодическую кусочно-линейную функцию  x + π + 2πk, 2 f (x) = −x + π + 2πk, 2 2πk − π ≤ x ≤ 2πk k ∈ Z. (39) 2πk ≤ x ≤ 2πk + π Будем предполагать, что в начальный момент времени значение фазовой переменной x1 лежит в области I : 0 < x1 ≤ π, а с течением времени переходит в область II : −π < x1 ≤ 0, в которой и заканчивает свое движение. Такое предположение обусловлено характеристиками модели и зависит от кон27

кретных ее параметров. Управление строится в виде кусочно-постоянно функции. Разобьем интервал времени T на N равных частей: t∗ = N ν, uj = u(t), t ∈ [(j − 1)ν, jν), j = 1, . . . , N . Учитывая вышесказанное и аппроксимацию (39), получим следующую задачу Zt∗ B(z) =C(z, t01 , t∗ ) = min  |u(t)| dt = min  u,t1 u,N1 0 N1 X |uj | + j=1  N X |uj | j=N1 +1 ẋ1 = A1 x1 + bu1 + c, x1 (0) = z, x11 (t1 ) = 0, t ∈ [0, t1 ], ẋ2 = A2 x2 + bu2 + c, x2 (t∗ ) = 0, t ∈ [t1 , t∗ ], (40) |u(t)| ≤ L, t ∈ T, где t01 — оптимальный момент перехода траектории x1 из области I в область II; матрицы A1 , A2 и векторы b и c вычисляются по формулам     θ1 1 0 0 0 θ1 1 0 0 0     θ2 − θ3 0 1 0 0 θ2 + θ3 0 1 0 0         1 2 A = θ4 − θ5 0 0 1 0 , A = θ4 + θ5 0 0 1 0 ,     θ − θ 0 0 0 1  θ + θ 0 0 0 1 7 7  6   6  −θ8 0 0 0 0 θ8 0 0 0 0 T πθ πθ πθ πθ 5 7 8 3 , , , . b = (0, 0, 0, 0, b0 )T , c = 0, 2 2 2 2 Воспользуемся формулой Коши общего решения для системы (40) Z N1 ν x1 (t1 ) =F 1 (N1 ν)z + F 1 (N1 ν − t)c dt+ (41) 0 + N1 X Z ∗ 1 1 F (N1 ν − t)b dt = k + uj j=1 2 jν (j−1)ν 2 N1 X uj dj , j=1 Z 1 Nν x (t ) =F (N ν − N1 ν)x (t1 ) + F 2 (N ν − t)c dt+ (42) N1 ν + N X j=N1 +1 Z jν 2 2 F (N ν − t)b dt = k + uj (j−1)ν N X j=N1 +1 28 uj dj ,

здесь k 1 и k 2 — константы, соответствующие первым двум слагаемым, для обоих равенств соответственно, dj , j = 1, . . . , N — множители при uj , а F 1 (t) и F 2 (t), t ≥ 0 — фундаментальные матрицы соответствующих однородных систем ẋ1 = A1 x1 и ẋ2 = A2 x2     1 1 F11 . . . F15 F11    .. . . . ..  =  ..  F 1 (t) =  .   .   . , 1 1 F51 . . . F55 F51  2 F11 2 F15  ...  .. . . . ..  F 2 (t) =  .   . . 2 2 F51 . . . F55 Перепишем задачу (40) с учетом решений (41) и (42) B(z) = min N X u,N1 x11 (t1 ) = |uj | j=1 F11 (N1 ν)z Z + N1 ν F11 (N1 ν − t)c dt+ (43) 0 + N1 X Z jν F11 (N1 ν uj j=1 x2 (t∗ ) = k 2 + − t)b dt = (j−1)ν N X k11 + N1 X uj d1j = 0, j=1 uj dj = 0, j=N1 +1 |u(t)| ≤ L, t ∈ T. Задача (43) — задача оптимального управления, в ней все элементы известны. Будем искать ее решение при помощи адаптивного метода Р. Габасова. Данная постановка задачи позволяет искать управления типа обратной связи, которое будем обозначать u(x(τ )), τ ∈ T , для оптимального управления и соответствующей ему оптимальной траектории введем обозначения u0 и x0 соответственно [15, 16]. Для того, чтобы в момент времени τ вычислить оптимальное управление u0ν (x0 (τ )), нужно решить задачу с начальными значениями z = x0 (τ ), с учетом того, что в момент времени τ −ν уже была решена задача с начальными значениями z = x0 (τ − ν). Кроме того на каждом шаге необходимо вычислять оптимальный момент перехода t01 , его вычисление 29

можно производить простым перебором, начиная со значения t01 (τ − ν). Если время на вычисление оптимального управления u0ν (x0 (τ )) и на оптимизацию параметра t1 не превосходит ν, то можно говорить о построении оптимального управления в режиме реального времени [17, 18]. Для фиксированных значений параметров модели построены оптимальное управление и оптимальные траектории. Результаты работы программы представлены на рис. 8–9. В данном примере оптимальным временем перехода из области I в область II является момент времени t01 = 0.266. 10 5 0 −5 −10 0 0.2 0.4 0.6 0.8 1 1.2 Рис. 8. Построенное управление 4 3 2 1 0 −1 −2 −3 −4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 Рис. 9. Интегральные кривые системы замкнутой построенным управлением 30

Заключение В данной работе рассмотрен метод построения оптимального управления, основанный на построении управления в виде кусочно-постоянной функции, поиск которой осуществляется после сведения задачи оптимального управления к задаче линейного программирования. Показано как можно адаптировать данный подход для разностных систем, неоднородных систем линейных дифференциальных уравнений и даже для некоторых классов нелинейных систем. Написан пакет прикладных программ в системе MATLAB, который реализует данный алгоритм. Программное обеспечение протестировано на примерах: рассмотрена задача стабилизации трехмассовой колебательной системы с оптимальным расходом топлива, а также аналогичная задача для двухмассовой системы с возмущениями; основываясь на реальных данных построена разностная динамическая модель межотраслевого баланса, проведено сравнение ее с дифференциальной моделью, построен прогноз, найдено оптимальное распределение инвестиций в секторы экономики; для нелинейной модели гибкого робота-манипулятора определено оптимальное управление. Результаты работы программы представлены на графиках. Стоит отметить, что, при проведении численных экспериментов, программный комплекс легко адаптировался к разным типам задач, а время работы программы было незначительным, учитывая реальные условия функционирования объектов управления. Также при построении оптимального управления для нелинейных моделей, а как следствие и для линейных моделей, возможно построение позиционного управления типа обратной связи. 31

Список литературы 1. Понтрягин Л. С., Болтянский В. Г., Гамкрелидзе Р. В., Мищенко Е. Ф. Математическая теория оптимальных процессов. М.: Наука, 1969. 384 с. 2. Беллман Р., Гликсберг И., Гросс О. Некоторые вопросы математической теории процессов управления. М.: ИЛ, 1962. 336 с. 3. Зубов В. И. Теория оптимального управления судном и другими подвижными объектами. Л.: Судпромгиз, 1966. 352 с. 4. Зубов В. И. Математические методы исследования систем автоматического регулирования. Л.: Машиностроение, 1974. 336 с. 5. Зубов В. И. Лекции по теории управения. М.: Наука, 1975. 496 с. 6. Калман Р. Е. Об общей теории систем управления. Труды I Междунар. конгресса ИФАК. М.: Изд-во АН СССР, 1961. Т. 2. С. 521–547. 7. Балашевич Н. В., Габасов Р., Кириллова Ф. М. Численные методы программной и позиционной оптимизации линейных систем управления// Журн. вычисл. математики и мат. физики. 2000. Вып. 40, № 6. С. 838–859. 8. Альсевич В. В., Габасов Р., Глушенков В. С. Оптимизация линейных экономических моделей. Минск: Изд-во БГУ, 2000. 211 c. 9. Габасов Р. Методы оптимизации: пособие. Минск: Четыре четверти, 2011. 472 с. 10. Федосеев В. В., Гармаш А. Н., Дайитбегов Д. М. и др. Экономикоматематические методы и прикладные модели: Учеб. пособие для вузов. М.: ЮНИТИ, 1999. 391 с. 11. Пересада В. П., Смирнов Н. В., Смирнова Т. Е. Управление развитием многопродуктовой экономики на основе динамической модели «затратывыпуск»// Вестник Санкт-Петербургского университета. Серия 10. При32

кладная математика, информатика, процессы управления, 2014. — № 4. — P. 121–134. 12. Попков А. С. Идентификация динамической модели межотраслевого баланса для экономики России и оптимальное распределение инвестиций на ее основе // Процессы управления и устойчивость. 2015. T. 2. № 1. С. 696– 701. 13. Краткие таблицы ресурсов и использования товаров и услуг [Электронный ресурс]: URL:http://www.gks.ru/free_doc/new_site/vvp/ tab-zatr-vip.htm (дата обращения: 13.03.15). 14. Krstic M., Kancllakopoulos I., Kokotovic P. Nonlinear and adaptive control design. New York: Library of Congress Catalog in Publication Data, 1995. p. 536. 15. Смирнов М.Н., Смирновa М.А. Вопросы синтеза стабилизирующих управлений при наличии неопределенных внешних возмущений // Процессы управления и устойчивость. 2015. T. 2. № 1. С. 503–508. 16. Смирнов М.Н. Метод учета ограниченных внешних воздействий при синтезе обратных связей с многоцелевой структурой // Вестник СанктПетербургского университета. Серия 10. Прикладная математика. Информатика. Процессы управления. 2014. № 2. С. 130–140. 17. Клюенков А. Л. Реализация адаптивного метода в одной задаче оптимального управления // Процессы управления и устойчивость. 2015. T. 2. № 1. С. 53–58. 18. Баранов О.В., Попков А.С., Смирнов Н.В. Оптимальная стабилизация квадрокоптера в режиме реального времени // Устойчивость и процессы управления: Материалы III международной конференции, посвященной 85-летию со дня рождения чл.-корр. РАН В.И. Зубова. СПб, 2015. С. 115– 116. 33

Приложение В приложении приведен код функций метода (MatLab). Вторая фаза адаптивного метода 1 f u n c t i o n [ J_bs , J_n , x , u , d e l t a , A1 , r ] = adapt_oc ( x , c , A, A1 , J_bs , J_n , d0 , d1 , u , d e l t a ) 2 [ n , N] = s i z e (A ) ; 3 [ nbs , q ] = s i z e (J_n ) ; 4 J_bs = s o r t ( J_bs ) ; 5 J = NaN(N, 1 ) ; 6 f o r j = 1 :N J( j ) = j ; 7 8 9 10 end %% локальные переменные r = 0; 11 l = z e r o s (N, 1 ) ; 12 t h e t a = NaN(N, 1 ) ; 13 14 15 J0 = NaN(N, 1 ) ; %% векторы оценок и потенциалов i f (N - nbs == n ) 16 A1 = i n v (A ( : , J_bs ) ) ; 17 u = c ( J_bs ) * A1 ; 18 d e l t a = NaN(N, 1 ) ; d e l t a ( J_n) = c ( J_n) - u *A ( : , J_n ) ; 19 20 21 22 end %% поиск оптимального направления и длины шага f o r j=J_n' i f ( d e l t a ( j ) < 0 ) && ( x ( j ) ~= d0 ( j ) ) 23 J0 ( j ) = d e l t a ( j ) ; 24 e l s e i f ( d e l t a ( j ) > 0 ) && ( x ( j ) ~= d1 ( j ) ) 25 J0 ( j ) = d e l t a ( j ) ; 26 e l s e i f ( d e l t a ( j ) == 0 ) && ( ( x ( j ) ~= d0 ( j ) ) 27 end 29 30 end 31 i f i s n a n ( J0 ) == 1 r = 1; 32 33 else 34 [ deltamax , j 0 ] = max( abs ( J0 ) ) ; 35 l ( j0 ) = sign ( delta ( j0 ) ) ; 36 l ( J_bs ) = ( - 1 ) * A1*A ( : , j 0 ) * s i g n ( d e l t a ( j 0 ) ) ; 37 f o r j = J_bs' if 38 l(j) > 0 t h e t a ( j ) = ( d1 ( j ) - x ( j ) ) / l ( j ) ; 39 elseif 40 l(j) < 0 t h e t a ( j ) = ( d0 ( j ) - x ( j ) ) / l ( j ) ; 41 else 42 theta ( j ) = i n f ; 43 end 44 45 end 46 if 47 | | ( x ( j ) ~= d1 ( j ) ) ) J0 ( j ) = d e l t a ( j ) ; 28 l ( j0 ) > 0 t h e t a ( j 0 ) = ( d1 ( j 0 ) - x ( j 0 ) ) / l ( j 0 ) ; 34

elseif 48 l ( j0 ) < 0 t h e t a ( j 0 ) = ( d0 ( j 0 ) - x ( j 0 ) ) / l ( j 0 ) ; 49 else 50 theta ( j0 ) = i n f ; 51 52 end 53 [ theta_s , j s ] = min ( t h e t a ( J_bs ) ) ; 54 [ theta_0 , j_0 ] = min ( [ t h e t a ( J_bs ( j s ) ) , t h e t a ( j 0 ) ] ) ; 55 x = x + theta_0 * l ' ; 56 i f j_0 == 2 % Замена опоры 57 b u f = s i z e (J_n ) ; 58 f o r j = 1 : buf ( 1 , 1 ) i f J_n( j ) == j 0 59 J_n ( j ) = [ ] ; 60 break ; 61 end 62 end 63 64 else 65 J_bs ( j s ) = [ ] ; 66 J_n = s e t d i f f ( J , J_bs ) ; 67 b u f = s i z e (J_n ) ; 68 f o r j = 1 : buf ( 1 , 1 ) i f J_n( j ) == j 0 69 J_n ( j ) = [ ] ; 70 break ; 71 end 72 73 end 74 J_bs = s e t d i f f ( J , J_n ) ; 75 end 76 end 77 end Первая фаза адаптивного метода 1 f u n c t i o n [ u , c , D, D1 , J_bs , J_n , d0 , d1 , v , d e l t a ] = t h e _ f i r s t _ p h a s e _ o c (D, b , c , d0 , d1 ) 2 [ n ,N] = s i z e (D ) ; 3 u = (D\b ) ' ; 4 5 w = b - D* u ' ; r = 0; 6 v = []; 7 i f w == z e r o s ( n , 1 ) d e l t a = [ ] ; D1 = [ ] ; 8 D = c a t ( 2 ,D, e y e ( n ) ) ; 9 d0 = c a t ( 2 , d0 , z e r o s ( 1 , n ) ) ; 10 d1 = c a t ( 2 , d1 , z e r o s ( 1 , n ) ) ; 11 c = cat (2 , c , zeros (1 ,n ) ) ; 12 u = cat (2 ,u , zeros (1 ,n ) ) ; 13 J = NaN(N+n , 1 ) ; 14 J_n = NaN(N, 1 ) ; 15 f o r j = 1 :N+n J( j ) = j ; 16 17 end 18 f o r j = 1 :N J_n ( j ) = j ; 19 20 end 21 J_bs = s e t d i f f ( J , J_n ) ; 22 w h i l e r == 0 35

[ J_bs , J_n , u , v , d e l t a , D1 , r ] = adapt_oc ( u , c , D, D1 , J_bs , J_n , d0 , d1 , v , d e l t a ) ; 23 24 [ nbs , q ] = s i z e (J_n ) ; 25 k = 0; 26 f o r j = 1 : nbs i f J_n ( j - k ) > N 27 J_n ( j - k ) = [ ] ; 28 k = k +1; 29 end 30 end 31 end 32 33 else 34 c1 = c a t ( 2 , z e r o s ( 1 ,N) , - 1 * o n e s ( 1 , n ) ) ; 35 D2 = e y e ( n ) ; 36 for i = 1:n i f w( i ) < 0 37 D2( i , : ) = ( - 1 ) * D2( i , : ) ; 38 end 39 40 end 41 D = c a t ( 2 ,D, D2 ) ; 42 u = c a t ( 2 , u , abs (w) ' ) ; 43 d0 = c a t ( 2 , d0 , z e r o s ( 1 , n ) ) ; 44 d1 = c a t ( 2 , d1 , abs (w) ' ) ; 45 46 J = NaN(N+n , 1 ) ; 47 J_n = NaN(N, 1 ) ; 48 f o r j = 1 :N+n J( j ) = j ; 49 50 end 51 f o r j = 1 :N J_n ( j ) = j ; 52 53 end 54 J_bs = s e t d i f f ( J , J_n ) ; 55 w h i l e r == 0 [ J_bs , J_n , u , v , d e l t a , D1 , r ] = adapt_oc ( u , c1 , D, D1 , J_bs , J_n , d0 , d1 , v , d e l t a ) ; 56 57 [ nbs , q1 ] = s i z e ( J_n ) ; 58 [ bss , q2 ] = s i z e ( J_bs ) ; 59 [ q3 , nn ] = s i z e ( u ) ; 60 i f u (N+1:nn ) == z e r o s ( 1 , nn -N) 61 c1 = c a t ( 2 , c , z e r o s ( 1 , nn -N ) ) ; 62 f o r j = 1: bss i f J_bs ( j ) > N && u ( J_bs ( j ) ) == 0 63 d0 ( J_bs ( j ) ) = 0 ; 64 d1 ( J_bs ( j ) ) = 0 ; 65 end 66 end 67 68 end 69 f o r j = 1 : nbs 70 if j > nbs break ; 71 72 end 73 i f J_n( j ) > N && u (J_n ( j ) ) < 1 e - 7 74 D = c a t ( 2 ,D ( : , 1 : J_n ( j ) - 1 ) ,D( : , J_n ( j )+1: nn ) ) ; 75 u = c a t ( 2 , u ( 1 : J_n( j ) - 1 ) , u ( J_n ( j )+1: nn ) ) ; 76 c1 = c a t ( 2 , c1 ( 1 : J_n ( j ) - 1 ) , c1 ( J_n ( j )+1: nn ) ) ; 77 d0 = c a t ( 2 , d0 ( 1 : J_n( j ) - 1 ) , d0 (J_n ( j )+1: nn ) ) ; 36

78 d1 = c a t ( 2 , d1 ( 1 : J_n( j ) - 1 ) , d1 (J_n ( j )+1: nn ) ) ; 79 f o r i = 1: bss 80 i f J_bs ( i ) > N && J_bs ( i ) > J_n ( j ) J_bs ( i ) = J_bs ( i ) - 1 ; 81 end 82 83 end 84 J_n ( j ) = [ ] ; nbs = nbs - 1 ; 85 end 86 end 87 end 88 89 end 90 D = D( : , 1 :N) ; 91 d0 = d0 ( 1 :N ) ; 92 d1 = d1 ( 1 :N ) ; 93 c = c ( 1 :N ) ; 94 u = u ( 1 :N ) ; 95 J = J ( 1 :N ) ; 96 for i = 1:n J_bs ( i ) = J ( i ) ; 97 98 end 99 J_n = s e t d i f f ( J , J_bs ) ; 100 end 37

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв