CC BY
85
Вестник МГТУ, том 18, № 1, 2015 г.
стр. 7-11
УДК 656.61.052.7
С.А. Агарков, С.В. Пашенцев
Параметрическая идентификация обобщенной модели Номото с помощью аппарата вариационного исчисления
S.A. Agarkov, S.V. Pashentsev
Parametric identification of the Nomoto generalized model using the apparatus of variational calculus
Аннотация. Предложен новый подход к идентификации параметров обобщенной судовой модели Номото с использованием аппарата классического вариационного исчисления и метода наименьших квадратов.
Abstract. A new approach to the identification of parameters of the Nomoto generalized vessel model has been proposed. The apparatus of classical calculus and the method of least squares have been used.
Ключевые слова: обобщенная модель Номото, вариационное исчисление, метод наименьших квадратов Key words: Nomoto generalized model, calculus of variations, method of least squares
1. Введение
Параметрическая идентификация модели является сложной математической задачей (Эльсгольц, 1969; Моисеев, 1979), которая может быть решена с помощью современной вычислительной техники.
Проблема идентификации параметров обобщенной модели Номото посредством разложения решений в ряды Фурье рассмотрена в статье (Пашенцев, 2010). Другой подход к идентификации предложен нами в работе (Yudin et al., 2014), где решена задача идентификации простейшей модели Номото для циркуляции судна. В настоящей статье предлагается вариационное решение задачи идентификации для обобщенного уравнения Номото (Nomoto et al., 1957) с использованием данных, полученных в ходе наиболее информативного стандартного испытания при выполнении маневра "зигзаг".
Дифференциальное уравнение, рекомендованное 14-й Международной конференцией опытовых бассейнов для решения проблемы управляемости (Соболев, 1976), определяет криволинейное движение судна и имеет вид
d2ra dra , , 3 _ ^ s ^ dO
Т —-—+ Т--+ га + к ra ra + v ra — KR о + KRT —
p dt2 s dt 1112 0 0 3 dt
(1)
где га - угловая скорость поворота судна; 5 - угол кладки руля; параметры Тр, Т,, T3, Ks, vb v2 подлежат идентификации по результатам натурных испытаний.
В дальнейшем будем оперировать уравнениями 1-го порядка; введем обозначения: Е = dra / dt -угловое ускорение судна; К - курс судна. Используем простое дифференциальное соотношение га = dК/dt. Таким образом, для описания рассматриваемого движения вместо уравнения (1) имеем следующую систему дифференциальных уравнений 1-го порядка:
dE
dt
_ ю _ к1ю|га|_ v2co3 + КоО + КоТ3 га0) / Тр,
dt
dt
dra
dt
dK
------— ra.
dt
(2)
Такое представление задачи дает возможность решать ее в вычислительной среде MathCad. Минимизируем следующий функционал:
min{j[a(K - K3)2 + (ra - o3)2]dt} = min{JAdt},
(3)
т.е. потребуем от модели максимальной адекватности экспериментальным данным по углу поворота судна и угловой скорости вращения. Интеграл (3) используем на интервале (0, tf). В качестве весового коэффициента выберем множитель a = (1 / f)2, делающий слагаемые однородными и равнозначными. Уравнение Эйлера - Лагранжа для экстремали в этом случае выглядит так:
7
Агарков СА., Пашенцев С.В. Параметрическая идентификация...
dF d ( dF Л |_ dF d / dF
дХ dt 1 v дХ' I dK" dt l . da
2a(K - K3) - d (2(a - шэ)) = 0, dt
что дает в итоге
da da3 n
a(K - K )----+----= 0.
dt dt
(4)
Если учесть, что da / dt = d2K/ dt2, то получим дифференциальное уравнение 2-го порядка относительно угла поворота курса судна К на экстремали:
d2 K da3
—-— aK = aK +------=w(t).
dt1 dt
Это уравнение решается известным образом; его общее решение записывается в форме
K (t )= E (t) e d + E2 (t) e
-t ltf
где E-i(t) и E2(t) находятся методом вариации констант в виде интегралов:
tf tj-
E (t) = I y(t)e -ttf dt, E2 (t) = I ffffe 1 E dt.
о о
В нашем случае можно также получить экстремаль в виде дифференциального уравнения 1-го порядка относительно угловой скорости судна ю. Подставим в уравнение (4) значение производной угловой скорости из уравнения (2), продифференцируем получившееся уравнение по времени и в него вновь введем значение производной угловой скорости. Получим нелинейное дифференциальное уравнение 1-го порядка относительно угловой скорости поворота на экстремали:
-<х(ю - К3) - юэ +(-Ts—~ со - VjCoU - v2co3 + Ks5 + К5Т3 —) IT = 0. (5)
dt dt F
Затем можно предпринять попытку идентифицировать параметры нашей модели. Перепишем уравнение (5) как линейное уравнение относительно параметров модели:
(а(ю-АГэ) + с6э)Г + rx(B + v1co|co|-i-v,2co3 -К§Ъ-КьТъЪ = -со.
(6)
Для этого следует иметь шесть условий, поскольку модель содержит шесть констант. С учетом начального условия на левом конце интервала интегрирования [ю(0) = 0, ю'(0) = 0 и 5(0) = 0] можно получить только одно алгебраическое уравнение, связывающее искомые константы:
аКэ(0) - ёэ (0) + К5Т3 /Т = 0.
р dt
На правом конце интервала интегрирования t = tf должно выполняться естественное граничное
условие:
dF dF
= 2(a - a3) = 0.
Откуда получим
дХ' да
ra(tf) = ю%).
Следовательно, необходимо набрать еще четыре условия для замыкания задачи идентификации. Зададим условия в промежуточных точках маневрирования. Тем самым из всех возможных экстремалей многопараметрического семейства выберем ту единственную экстремаль, которая пройдет через заданные точки в пространстве состояний модели. Вычислительные процедуры опишем ниже в процессе численного решения конкретной задачи идентификации.
2. Численное решение
Рассмотрим численный пример использования этого подхода аналогично тому, как это было сделано при решении более простых задач разгона и циркуляции судна (Yudin et al., 2014). Результаты натурных испытаний типа "зигзаг", как правило, слишком "зашумлены" погрешностями различного генезиса. Поэтому в качестве опытных данных используем результаты моделирования маневра "зигзаг 10 /10" танкера "Саратов" (в балласте). Параметры математической модели указаны в таблице.
T 1 p T 1 s V: V2 Ks T3
291 11 -133 6 815 0.0285 0.114
8
Вестник МГТУ, том 18, № 1, 2015 г.
стр. 7-11
Данные параметры были получены при решении задачи идентификации обобщенного уравнения Номото посредством разложения уравнения движения в ряд Фурье (Пашенцев, 2010).
Результаты такого моделирования представлены на рис. 1: курс судна, угловая скорость поворота и кладка руля показаны как функции времени.
о- ■ X -ч у
\ :/~У “Л \ ■ ,д \ -
о- Ч'. —Ч.
60*-------------------------1-------------------------1-------------------------1-------------------------1-------------------------
0 100 200 300 400 500
t
время в сек
-----Кладка руля del
----- Курс К
----- Угл. скорость Ош
Рис. 1. Изменение характеристик движения судна при испытаниях "зигзаг 10 / 10"
В качестве дополнительных десяти точек возьмем характеристики движения, достигающие экстремальных или нулевых значений (рис. 1). Набор значений характеристик, полученных с помощью режима трассировки в системе MathCad, представлен на рис. 2. Вектор Т определяет моменты времени, в которые требуется совпадение характеристик состояния модели и результатов испытания; вектор В -правые части уравнения (6), матрица А содержит коэффициенты при искомых параметрах модели в том же уравнении. При этом система десяти линейных уравнений записывается в матричной форме:
AX = B,
где X - вектор искомых параметров.
' 5 ' 50 54
1
1 -1.693Т0*5
2 -8.20Г10-3
3 -9.002-10-3
4 0.015
5 0.015
6 -0.014
7 0.016
8 0.014
9 0.011
10 -0.016
1 2 3 4 5
8.385-10-6 1.266-10-5 2.867-10-М 4.854 Т0-15 -0.175
2 -1.132-10-6 2.14-10-1 6.725-10-5 5.515-10-7 -0.131
3 -1.385-10-5 1.835-10-4 8.103-10-5 7.294*10-7 0.044
4 1.373-10-5 -7.105-10-5 -2.29-10-4 -3.465*10-6 0.175
5 1.33-10-5 7.57Г10-5 -2.248-10-4 -3.37*10-6 0.044
6 -1.562-10-5 -1.773-10-4 2.055-10-4 2.945-10-6 0.044
7 1.737-10-5 -7.433-10-5 -2.47Г10-4 -3.884-10-6 0.175
8 2.736-10-5 2.643-10-4 -1.866-10-4 -2.548-10-6 -0.175
9 7.6Т0-6 -2.957*10-4 -1.288Т0-4 -1.463-10-6 0.175
10 -1.762-10-5 8.2Г 10-5 2.473-10-4 3.889-10-6
Рис. 2. Исходные данные,
используемые для определения параметров модели
Данные значения характеристик состояния объекта (угловая скорость, курс, угол кладки руля) выберем в качестве экспериментальных (они нами наблюдаемы). При этом учтем, что в условиях натурного эксперимента обычно не наблюдаются значения производных курса и угловой скорости, которые необходимы для подстановки в уравнения (5) и (6). Их найдем с помощью конечных разностей по традиционным формулам. Следует отметить, что на рис. 2 приведено избыточное количество данных для однозначного решения задачи. Это обстоятельство предоставляет возможность решить задачу определения параметров как переопределенную и получить большую устойчивость решения.
Дальнейшее решение можно осуществить с помощью MathCad двумя способами. Использование встроенной функции lsolve дает решение немедленно (рис. 3, скриншот экрана).
' 191.713 ' 4.98
-136.178
6.112 х 103 0.016
, -0.028 ,
Рис. 3. Решение переопределенной системы с помощью встроенной функции lsolve
9
Агарков СА., Пашенцев С.В. Параметрическая идентификация...
Другое решение требует определения для матрицы А псевдообратной матрицы Арр (рис. 4).
Д» JaT.A}-1.AT
.waiiRv '
Х^:= Арр ■ В
' 117393 " 1.409 -137.677
5.571 х 10J
S.48 х 10” 3 , -002 ,
( 291 > 11
-133 6815 0.02S5 0.114 j
Рис. 4. Решение переопределенной системы с помощью псевдообратной матрицы Арр
Решения Х и Х получились совершенно одинаковыми, значит, функция lsolve использует псевдообратную матрицу. На рис. 4 приведен также вектор параметров, которые послужили базой для получения опытных данных в модели Номото.
Наконец, можно максимально переопределить задачу и использовать весь комплекс измерений кинематических параметров с номерами от 1 до 500 (500 с - длительность эксперимента, рис. 1). В данном случае естественно применить метод наименьших квадратов (МНК) и получить матрицу А при искомых параметрах и вектор правой части В. Для нашей системы уравнений формальное применение МНК состоит в умножении уравнения (6) в точке с номером к последовательно на ю'', ю', ю|ю|, ю3, 5, 5', вычисленные также в точке с номером к, затем в сложении по всем экспериментальным точкам. Получим так называемую нормальную систему шести уравнений с шестью неизвестными, решение которой проводим обычным образом, демонстрируя его в виде фрагмента решения в среде MathCad (рис. 5). Чтобы отличить данное решение от предыдущих, введем обозначения: С - матрица системы; D - вектор свободных частей. Оба решения (прямое и полученное с помощью псевдообратной матрицы) практически совпадают. Как и ожидалось, результат ближе к параметрам модели, использованной для генерации экспериментальных данных.
Срр := (сТ с) • СТ X]^:= Срр D
Х2 := С~ 1 - D
' 240.345 " 8.177 -135.723
6.572 х 103 0.023
. -0.037 .
4'
( 291 " И
-133 6815 0.0285 I 0.114,
Рис. 5. Решение системы нормальных уравнений методом наименьших квадратов
3. Выводы
Предложен вариационный подход к решению задач параметрической идентификации математических моделей. При моделировании с учетом малого числа параметров такой подход способствует получению достаточно точных результатов (Yudin et al., 2014) и может быть использован для аппроксимации сложных движений набором простейших движений, модели которых легко идентифицируются. В настоящей статье такой подход был применен к модели с большим числом идентифицируемых параметров и дал вполне удовлетворительный по точности результат. Наилучший результат идентификации получен при использовании всех исходных данных модельного эксперимента, т.е. при двойном отборе идентифицируемых параметров - с помощью вариационного уравнения и метода наименьших квадратов.
Литература
Nomoto K., Taguchi T., Honda K., Hirano S. On steering qualities of ships. JSP. 1957. N 35. P. 56-64.
Yudin Yu., Pashentsev S., Petrov S. Using Pontryagin maximum principle for parametrical identification of ship maneuvering mathematical model. Transport Problems. 2014. V. 9, Issue 2. P. 11-18.
Моисеев Н.Н. Численные методы синтеза оптимальных управлений. М., Наука, 1979. C. 443.
Пашенцев С.В. Параметрическая идентификация маневренных характеристик по результатам испытаний типа "Зигзаг". Вестник МГТУ. 2010. Т. 13, № 4/1. С. 730-735.
Соболев Г.В. Управляемость корабля и автоматизация судовождения. Л., Судостроение, 1976. С. 478. Эльсгольц Л.Э. Дифференциальные уравнения и вариационное исчисление. М., Наука, 1969. С. 375.
10
Вестник МГТУ, том 18, № 1, 2015 г.
стр. 7-11
References
Nomoto K., Taguchi T., Honda K., Hirano S. On steering qualities of ships. JSP. 1957. N 35. P. 56-64.
Yudin Yu., Pashentsev S., Petrov S. Using Pontryagin maximum principle for parametrical identification of ship maneuvering mathematical model. Transport Problems. 2014. V. 9, Issue 2. P. 11-18.
Moiseev N.N. Chislennye metody sinteza optimalnyh upravleniy [Numerical methods for the synthesis of optimal controls]. M., Nauka, 1979. P. 443.
Pashentsev S.V. Parametricheskaya identifikatsiya manevrennyh harakteristik po rezultatam ispytaniy tipa "Zigzag" [Parametrical identification of maneuvering characteristics according to results of environmental tests of the "Zigzag" type in non-linear model of vessel controllability]. Vestnik MGTU. 2010. T. 13, N 4/1. P. 730-735.
Sobolev G.V. Upravlyaemost korablya i avtomatizatsiya sudovozhdeniya [Ship handling and automation of navigation]. L., Sudostroenie, 1976. P. 478.
Elsgolts L.E. Differentsialnye uravneniya i variatsionnoe ischislenie [Differential equations and the calculus of variations]. M., Nauka, 1969. P. 375.
Информация об авторах
Агарков Сергей Анатольевич - Институт экономики, управления и права МГТУ, кафедра управления социально-экономическими системами, д-р экон. наук, профессор
Agarkov S.A. - Institute of Economics, Management and Law MSTU, Department of Management of Social and Economic Systems, Dr of Econ. Sci., Professor
Пашенцев Сергей Владимирович - Морская академия МГТУ, кафедра судовождения, профессор, канд. физ.-мат. наук, e-mail: serpass15@yahoo.com
Pashentsev S.V. - MSTU Marine Academy, Navigation Department, Professor, Cand. of Phys.-Math. Sci., e-mail: serpass15@yahoo.com
11
