CC BY
44
Идентификация некоторых маневренных характеристик судна по результатам натурных испытаний
С.В. Пашенцев
Судоводительский факультет МА МГТУ, кафедра судовождения
Аннотация. Предлагается два подхода к идентификации ряда маневренных характеристик циркуляции судна (радиус, угловая скорость поворота) по результатам натурных испытаний с использованием метода наименьших квадратов (МНК). Один из них относится к обработке результатов, полученных на тихой воде, второй - к испытаниям в ветровых условиях. Последнее особенно важно, поскольку почти невозможно провести натурные испытания при идеально тихой погоде.
Abstract. There are two ways to identify the number of maneuvering characteristics of the vessel's circulation (radius, rate of turn) as a result of environmental tests using the least-squares method. One of them belongs to processing of the results obtained on still water, the second - to tests at windy conditions. The last is especially important, as it is impossible to have any environmental tests at ideally calm weather.
1. Введение
Проблема идентификации маневренных характеристик циркуляции не столь проста, как это кажется с первого взгляда. Хотя установившаяся циркуляция и входит в стандартные испытания судна, но определяемые по ее результатам характеристики, прежде всего радиус установившейся циркуляции, содержат существенные погрешности. Они связаны с обработкой "на глаз" результатов испытаний, при этом весь объем траекторной информации, полученной в ходе испытаний, практически не используется. Вторая проблема связана с невозможностью провести испытания в идеальных гидрометеорологических условиях, поскольку на любой акватории действует ветер той или иной силы. Именно эти проблемы и решает соответствующая обработка результатов для получения маневренных характеристик судна на циркуляции, которая использует весь объем траекторных измерений, проведенных в ходе испытаний. Это позволяет научно обоснованным способом находить важнейшие маневренные характеристики судна.
2. Идентификация радиуса и центра установившейся циркуляции в отсутствие ветра
Пусть проведены испытания судна на циркуляции, первично обработаны результаты в виде данных GPS (лучше DGPS!) и найдены координаты точек траектории, полученных судовой системой обсервации. Данные переведены из навигационных координат в обычные декартовы координаты и представлены в метрах. Такая первичная обработка не представляет никаких сложностей и легко реализуется на любом компьютере. Особенно просто делается это, если система GPS соединена напрямую с компьютером; тогда данные обрабатываются в реальном времени.
Предположим, что мы имеем координаты траектории как пары (Xk, Yk) в количестве N и нам следует получить на их основании радиус циркуляции Rc и координаты ее центра (Xc, Yc). Для каждой точки траектории (Xk, Yk) можно найти квадрат расстояния до окружности циркуляции с уравнением (X - Xc)2 + (Y - Yc)2 = Rс2
dk = (Xk - Xc)2 + (Yk - Yc)2 - Rc2, (1)
причем параметры окружности нам неизвестны и их следует определить по траекторным измерениям.
Естественно, что точки, полученные в испытаниях, не попадают на окружность за счет погрешностей измерений и других источников не идеального проведения реальных испытаний (течение, ветер, волнение, рыскание судна). Поэтому сформулируем задачу выбора параметров окружности (Xc, Yc) и Rc как задачу минимизации суммы квадратов отклонений точек реальной траектории от выбранной окружности:
min{F(Xc, Yc, Rc)} = Min{Z dk2} = min{Z [(Xk - Xc)2 + (Yk - Yc)2 - Rc2 ]2}, (2)
где сумма берется по всем точкам или выбранной части точек натурной траектории, а минимум ищется вариацией искомых параметров циркуляции Xc, Yc и Rc. Процедура минимизации формально проста и требует отыскивать точку минимума функции трех переменных, приравняв к нулю три частных производных:
dF/dXc = 0, dF/dYc = 0, dF/dRc = 0. (3)
Реализация частного дифференцирования и освобождение от несущественных при равенстве нулю постоянных множителей приводит к трем нормальным уравнениям, характерным для подобных проблем минимизации:
2 [(X - X) + (Ук - Ус)2 - Яс2] = 0,
2 Хк [(Хк - Хс)2 + (Ук - Ус)2 - Яс2] = 0, (4)
2 Ук [(Хк - Хс)2 + (Ук - Ус)2 - Яс2] = 0.
Первое из этих уравнений легко приводит к соотношению:
Яс2 = 2 [(Хк - Хс)2 + (Ук - Ус)2 / N (5)
которое можно подставить во второе и третье уравнения и получить
2 (Хк - Хср) [(Хк - Хс)2 + (Ук - Ус)2] = 0, 2 (Ук - Уср) [(Хк - Хс)2 + (Ук - Ус)2] = 0,
где введены средние арифметические значения обработанных координат Х и У:
Хср = (2 Хк) / N Уср = (2 Ук) / N. (7)
Фактически (6) уже представляют два уравнения относительно двух неизвестных координат центра окружности, но их явное решение не слишком удобно, и потому проведем предварительно некоторые преобразования. Запишем эти уравнения только в отклонениях координат от средних значений. Это упрощает сами выражения и сделает ниже порядок входящих в них выражений, что важно при численной реализации задачи. Отклонения запишем как хк и Ук, не вводя новых имен переменных. В терминах отклонений уравнения (6) получат вид:
2 (Хк) [(Хк - Хс)2 + (Ук - Ус)2] = 0, (8)
2 (Ук) [(Хк - Хс)2 + (Ук - Ус)2] = 0,
где под Хс и ус также понимаются отклонения координат центра от средних значений Х и У. Комбинируя эти два уравнения, можно выразить Ус через Хс линейно, хотя сами уравнения (8) относительно этих параметров квадратичны:
Ус = АХс + В, (9)
где коэффициенты А и В выражаются через ряд сумм, вычисленных по отклонениям координат. Вначале введем эти суммы, обозначая их именными индексами суммируемых координат:
„ _ ¿хх ' = 2Хк 'Хк _ 3 „ _ 2 „,
¿х 2Хк? V"1 ^ххх 2Хк > ¿хху 2Хк 'Ук /"1 ПЛ
5ху = 2ХкУк, ^ 3 V 2 (10)
^у = 2Ук, / = ~ , У 5ууу = 2Ук , ^хуу = 2ХкУк .
вуу = АУк'Ук,
Эти суммы позволяют выразить коэффициенты А и В так:
А = (5ху-5х - ¿хх'^у) / (5ху-5у - 5уу-5х)
В °.5[(*ххх + ^хуу)^у (^ууу + 5хху)5х)] / (5ху-5у ¿уу'^хХ
(11)
Подставляя Ус в виде (9) в любое из двух уравнений, получим единственное квадратное уравнение относительно координаты Хс. Ниже даем результат подстановки в первое из уравнений (8) (для проверки это было сделано и для второго уравнения с получением совпадающего результата):
Хс25х(1 + А2) - 2Хс(5хх + А^ху - А-В-^х) + (5ххх + ^хуу - 2В-5ху + В^х) = 0. (12)
Далее решение очевидно - по найденному значению хс находим Ус по (9) и, наконец, Яс по (5).
3. Численная обработка данных модельного и натурного экспериментов
Несмотря на принципиальную несложность формул (5-12), реализовать их эффективно, особенно в судовых условиях, можно только с помощью вычислительной техники. Первой проверкой правильной работы этих формул явилась проверка модельная, когда генерируется идеальная траектория в виде окружности, координаты точек на ней "портятся" погрешностями некоторого уровня с и затем по этим искаженным точкам восстанавливаются координаты центра окружности и ее радиус. Для этого была разработана специальная исследовательская программа (Пашенцев и др., 2003), включавшая в себя, в частности, описанную процедуру идентификации. Результаты работы программы приведены в табл. 1. В ней первая колонка содержит число сгенерированных точек, равномерно расположенных по окружности. Вторая колонка содержит величину средней квадратической погрешности сгенерированных точек при их нормальном распределении. Третья, четвертая и пятая колонки содержат координаты центра окружности и ее радиус, которые были определены по формулам (5-12). Первая строка тела таблицы дает значения этих параметров, которые были использованы при генерации и с которыми мы
сравниваем вычисляемые значения. Видно, что точность восстановления параметров по координатам траектории зависит как от внесенной погрешности точек, так и от числа точек. Чем больше погрешность сгенерированных точек и чем меньше самих точек, тем больше отклонения найденных значений параметров от их истинных величин, что и демонстрирует табл. 1.
К сожалению, имея на руках большое число результатов натурных экспериментов, мы не можем использовать их в данной методике обработки с гарантией получения надежных выводов. Дело в том, что все эксперименты проводились в условиях ветра, и потому движение судна проходило по кривым, далеким от окружности. Кроме того, результаты обработки фактически не с чем сравнивать, в первую очередь это касается координат центра циркуляции. Однако все же приведем результат такой обработки одной из экспериментальных циркуляций, проведенной капитаном Кислым на судне "Шаньков" АО Тралфлот в 2003 г. согласно рекомендациям (Кацман и др., 1970). Результаты такой обработки представлены на рис. 1, где условия эксперимента и полученные значения параметров циркуляции приведены в подрисуночных надписях. Из рисунка хорошо видно, что результат существенно зависит от того, какие точки из числа экспериментальных будут обрабатываться. Это же в значительной мере субъективно, потому что трудно отличить, например, эволюционный период циркуляции от ветрового сноса, что и имеет место в данном примере. Именно поэтому ниже приведена методика идентификации, разработанная автором для учета действия ветрового дрейфа.
Таблица 1
N о, м Хс, м Гс, м Яс, м
0 20 -30 100
30 10 25.12 -41.09 102.85
60 10 22.68 -35.89 101.21
240 10 20.05 -30.47 100.79
60 5 22.25 -32.67 100.25
360 15 19.14 -30.24 103.14
Рис. 1. Идентификация параметров натурной циркуляции (ветер 10 м/с 90 град., руль 15 град. п/б) а) - по всему массиву точек эксперимента (58 точек, пунктирная окружность: Хс = -78 м,
Гс = -54 м, Яс = 52 м); б) - по урезанному массиву точек эксперимента (только закрашенные точки с 5 по 43, сплошная окружность: Хс = -69 м, Ус = -61 м, Яс = 44 м)
4. Идентификация радиуса установившейся циркуляции и угловой скорости поворота в условиях действия ветра
Из сказанного выше следует, что идентификация параметров циркуляции в реальных условиях ветрового воздействия есть задача более сложная. Однако ее можно решить, также применяя метод наименьших квадратов, но не к круговой траектории, а к некоторой спиралевидной траектории. Именно по такой траектории движется судно, когда оно выполняет циркуляцию при одновременном ветровом дрейфе. Опишем такую траекторию, как одновременное движение по окружности, при этом сама окружность сносится с некоторой скоростью Ук в некотором направлении а. Можно оперировать также двумя проекциями скорости Ух и Уу вместо самой скорости и ее направления. При этом подходе мы будем иметь четыре требующих определения параметра - радиус циркуляции Яс, угловая скорость циркуляции а, компоненты скорости ветрового дрейфа Ух и Уу.
Введены одновременно параметры Яс и а, но при циркуляции они связаны простым соотношением а -Яс = У, где У - линейная скорость на циркуляции. Но на циркуляции даже при неизменной работе двигателя судна линейная скорость падает и потому не может считаться постоянной.
Поэтому более гибким будет подход независимого варьирования угловой скорости и радиуса циркуляции. После их определения можно будет определить еще и величину линейной скорости при повороте.
Подобно (2), потребуем минимума квадрата отклонений опытных точек траектории судна (Хк, Ук) от траектории сносимого ветром круга циркуляции:
шВДУх, V,, Яс)} = ш1п{2 [(Яс - Яс сов(оТ*) + УхТк -Хк)2 + (Яс в1п(оТ*) + УуТк - Г*)2]}. (13)
Как и ранее, потребуем равенства нулю четырех частных производных
дПдЯс = 0, дПда = 0, дПд¥х = 0, дПд¥у = 0, (14)
что приводит к четырем нормальным уравнениям МНК.
Развернутые выражения для производных в том же порядке таковы:
2[(1 - сов(оТк))(Яс - Яс сов(оТ*) + УхТк - Хк) + ьЩаТк) (Яс ¡ип(юТк) + УуТк - Г*)] = 0, 2 [(Тк ып(аТк))(Яс - Яс сов(аТ*) + УХТ* - Хк) + Т* сов(аТ*) (Яс в1п(аТ*) + УуТ* - Г*)] = 0, (15)
2 [Тк(Яс - Яс сов(®Т*) + Ух Тк - Хк)] = 0, 2 [Тк(Яс ът(тТк) + У, Т* - Г*)] = 0.
Представим производные в форме, содержащей ряд сумм, которые обозначим, следуя принципу, введенному выше. Так, "сов означает сумму 2Тк• сов(о -Тк), ^4х81п означает 2ТкХк- вш(о • Тк) и т.п. В этих обозначениях система уравнений (14) будет иметь следующий вид:
(15а)
(16)
Яс(Х - 2"сов) + Ух(" - + Уу• "[зт + "хсов - "уВ1п = 0, Яс("[ - "сов) + Ух• - "х = 0,
Яс'"вт + Ух• "шип + Уу• "[[сов - "[хвт - "усов = 0, Яс'"вт + Уу• - "у = 0.
Легко выразить проекции скорости ветрового сноса из третьего и четвертого уравнений (15а):
Ух = ["х - Яс• (" - "сов)] / "ь Уу = ("у - Яс • "1з1п) / "б
и подставив их в уравнения 1 и 2 той же системы, получим систему для Яс и о :
Яс{2(Х - "сов!" - (" - ";сов) - О^ЬшО } = {(-"хсов + "увт^й - "х(" - - "у•"втК (17)
Яс{"зт- (" - "сов^шт - "вШ"[сов} = {("хвт + "усов^и - "х • "зт - "V"[сов}.
Выразив Яс из каждого уравнения (17) и приравняв результаты, получаем, наконец, уравнение, которое содержит только один искомый параметр - угловую скорость циркуляции о :
{(-"хсов + "увт^й - "х (" - "сов) - "V"зт} / {2(Х - "сов^й - (" - "сов) - О^з^ } = (18)
= {("хв1п + "усов)": - "х• "ивш - "[у•"сов} / {"зт- (" - "сов^Шт - "¡т•"сов}.
Запишем и выражение для Яс, полученное из первого уравнения (17), т.к. оно послужит для нахождения Яс после того, как будет определена угловая скорость:
Яс = {(-"хсов + "вт!" - "х(" - "[сов) - "у•"вт} / {2(Х - "сов!" - (" - "[сов) - ("цг^ } = / Н\. (19)
Аналогично находится Яс из второго уравнения (17):
Яс = {("хвт + "усов!" - "х • "[^¡п - "[у• "[сов} / - (" - "сов^йвт - "зт•"сов} = ^2 / Н2. (19а)
Уравнение (18) есть сложное трансцендентное уравнение, которое можно решать только численно, используя метод последовательных приближений. Это связано с тем, что входящие в (18) суммы можно найти, только предварительно задав угловую скорость. Впрочем, метод дихотомии работает здесь безупречно.
Хотя мы рассмотрели случай независимых величин Яс и о, все же проделаем процедуру минимизации для более естественного случая, когда эти характеристики связаны простым соотношением. В этом случае поиск минимума из безусловного переходит в условный с одним дополнительным условием в форме о Яс - У = 0. Процедура поиска минимума немного усложняется - вводится коэффициент Лагранжа Я и ищется безусловный минимум для модифицированной следующим образом функции:
ш1п{^(Ух, Уу, Яс) - Я(о Яс - У)}, (20)
вариацией тех же переменных. При частном дифференцировании по этим четырем переменным мы получим систему (15) с небольшими изменениями. Третье и четвертое уравнения системы вообще не изменятся, т.к. введенное условие не зависит от переменных Ух и Уу. Первое и второе уравнения
изменятся минимально: первое уравнение системы (15) будет иметь в правой части член Х-ю, а второе -член XRc. Это замечание относится и к другой форме записи системы (15а). Отсюда следует, что выражения (16) для Vx и Vy останутся неизменными, а два первых уравнения получат форму:
Rc(N - 2Scos) + Vx(St - Scos) + Vy Stsin + Sxcos — Sysin = Хю, (21)
Rc-Stsin + Vx'Sttsin + Vy'Sttcos - Stxsin — Stycos = XRc.
Каждое из этих уравнений определяет коэффициент Лагранжа. Приравняв оба выражения для Х и учтя, что Rc = V/ю, вновь получим одно трансцендентное уравнение для неизвестной угловой скорости ю. Его можно записать, используя уже известные выражения, введенные в (19) и (19а):
®2(V • H2 - ю -G2) = V( V • H1 - ю -Gj). (22)
Решаем его численно, после чего легко находится Rc = V/ю, а с помощью соотношений (16) -компоненты скорости ветрового дрейфа Vx и Vy.
Еще один важный вариант минимизации - при двух дополнительных концевых условиях равенства координат модельной и идентифицирующей кривых. Это типичная задача минимизации с краевыми условиями или, иначе, двухточечная краевая задача. Условия эти выглядят так:
Rc - Rccos(^TN) + Vx-TN = XN, (23)
Rc sin(<® -Tn) + Vy-TN = Yn,
где TN означает момент времени, в который уравниваются координаты. При этом в левых частях равенств (23) расположены координаты идентифицирующей кривой, а в правой части - модельной. Выберем момент уравнивания координат, соответствующий полному циклу движения по окружности, т.е. TN = 2п/ю. В этом случае выражения (23) существенно упрощаются и получают форму
Vx-Tn - Xn = 0, Vy-TN - Yn = 0, (24)
которую мы и используем при конструировании новой функции для минимизации. Она содержит два коэффициента Лагранжа соответственно двум дополнительным условиям (24):
min{F(Vx, Vy, Rc) - Xx(Vx-Tn - Xn) - X2(VX-TN - XN)}. (25)
Теперь при выполнении частного дифференцирования два первых уравнения системы (15, 15а) не изменятся, т.к. дополнительные условия не содержат переменных Rc и ю. Два же следующих уравнения, которые возникают при дифференцировании по Vx и Vy, немного видоизменятся, получив в правых частях члены Xx-TN и X-TN вместо нулей:
Rc(St - Stcos) + Vx'Stt - Stx = Х1 • TN (26)
Rc-Stsin + Vy-Stt - Sty = X2-TN.
Эти два уравнения позволяют найти коэффициенты Лагранжа, которые необходимы для вычисления самого минимума минимизируемой функции. Но нас интересуют другие параметры, и потому из условий (24) найдем Vx и Vy и подставим их в два первых уравнения (15а), которые не изменились. Получим
Rc(N - 2Scos) + XNiSt - Scos) / TN + YN -Stsin / TN + Sxcos - Sysin = ° (27)
Rc'Stsm + XN 'Sttsm / TN + YN -Sttcos / TN - Stxsin - Stycos = 0.
Выразив из каждого уравнения (27) радиус Rc и приравняв эти выражения, получим новое трансцендентное уравнение для угловой скорости ю :
[XN(St - Scos) / TN + YN '^tsin / TN + Sxcos - Sysin] / (N - 2Scos) = (2g)
= (XN -Sttsin / TN + YN 'Sttcos / TN - Stxsin - Stycos) / Stsin.
Решив его численно, находим затем Rc из любого из уравнений (27). Компоненты скорости уже известны, и тем самым решение задачи завершается. Численный анализ показывает, что этот вариант минимизации дает идентификационную траекторию, наиболее близкую к модельной, но при этом сами параметры идентифицируются хуже. Это и определяет сферу применимости каждого метода идентификации из предложенных выше.
5. Численная идентификация параметров в условиях действия ветра
Как и в разделе 4, применим разработанные методы идентификации вначале к модельным результатам, а затем и натурным. В отличие от случая чистой циркуляции, генерировать модельное
движение будем не чисто кинематически, а используем имеющуюся в нашем распоряжении динамическую модель (Юдин, Пашенцев, 2002) судна (танкер "Саратов"), разработанную автором совместно с Юдиным Ю.И. в рамках хоздоговорной работы с АО ММП в 2003 г.
Моделировалось движение танкера в условиях действия ветра скорости 16 м/с и направления 60 град., при этом судно совершало циркуляцию с начальным направлением на N со скоростью 4 уз. На рис. 2 представлены три кривые, которые характеризуют результаты идентификации. Кривая с номером 1 описывает циркуляцию модели в отсутствие ветра с кладкой руля 35 град. п/б. Эта циркуляция имеет радиус Яс = 279 м и угловую скорость 24.66 град./мин., что хорошо согласуется с результатами сдаточных испытаний танкера. Они проводились на Балтике летом 2002 г., ветер практически отсутствовал. Кривая 3 описывает движение модели в условиях ветра, о котором сказано выше. Кривая 2 дает идентификационную траекторию движения, при этом получены такие параметры идентификации: радиус Яс = 271 м, угловая скорость о = 24.4 град/мин. Это соответствует скорости на циркуляции V = 1.923 м/с вместо 2.06 м/с, которые следуют из исходной скорости судна в 4 уз., что соответствует известному факту падения скорости на циркуляции. Видим, что получены близкие результаты по двум основным параметрам. Кроме того, получена и дополнительная информация - судно сносится ветром в направлении 198 град. со скоростью 0.47 м/с. Заметим, что кривые 2 и 3 рис. 2а не совпадают локально, но поскольку глобальные параметры совпадают, можно считать подобный подход вполне оправданным. Ведь нашей целью было определение именно этих глобальных параметров. Однако это несовпадение траекторий наводит на мысль дальнейшего уточнения кинематической модели, скажем, снос окружности с одновременным изменением радиуса. Это есть основание для продолжения исследований по идентификации в условиях действия ветра.
Численные идентификации были проведены для различных направлений действующего ветра. В зависимости от этого фактора и параметры идентифицируются различным образом. Результаты расчетов приведены в табл. 2 в зависимости от направления ветра Ка: радиус Яс, угловая скорость о, и параметры дрейфа - скорость и направление. Видно, что наилучшие, с нашей точки зрения, результаты идентификации при направлении ветра порядка 60 град. Именно такой курсовой угол ветра можно рекомендовать в натурных экспериментах с целью последующей обработки предложенным методом для получения радиуса "чистой" циркуляции. Неплохой результат дает также осреднение радиуса, полученного для разных курсовых углов ветра (Яс)ср = 260 м. Попутно полученные параметры дрейфа судоводитель может использовать для выполнения маневра при действии ветра.
Для нулевого направления ветра, когда судно начинает разворачиваться точно на ветер, результаты в виде кривой идентификации приведены на рис. 2б, который можно качественно сравнить с рис. 2а. В обоих случаях отметим еще раз локальное несовпадение кривых модельной и идентифицирующей. Для случаев, когда важна большая степень локального соответствия, использована для расчета методика совпадения конечных точек траекторий. Этот случай приводится на рис. 2в при ветре направления 60 град. При этом получены несколько иные результаты идентификации: Яс = 285 м, о = 24 град./мин., Vдр = 0.47 м/с, Кдр = 207 град., что достаточно хорошо совпадает с результатами табл. 2.
Таблица 2
Ка, град. 0 30 60 90 120 150 180
Я„ м 290 286 271 255 245 239 236
о, град./мин 22.41 23.05 24.39 26.14 27.85 29.25 29.85
Vдр. м/с 0.33 0.39 0.47 0.48 0.44 0.35 0.23
Кдр. град. 144 176 198 216 235 258 290
Рис. 2. Идентификация параметров модельной циркуляции при ветре (кривые: 1 - циркуляция без ветра, 2 - идентифицирующая, 3 - модельная)
Пашенцев С.В. Идентификация некоторых маневренных характеристик судна... 5. Заключение
Предложенные методы определения маневренных характеристик судна по результатам натурных испытаний в отсутствие ветра и при его действии показали, что они обладают необходимыми качествами для адекватного определения радиуса и угловой скорости циркуляции. Знание именно этих характеристик важно при непосредственном выполнении маневра. Показано, что различная постановка задачи минимизации функционала влияет на определение указанных характеристик примерно в той мере, которая определена для допустимой погрешности их нахождения - порядка 10 %. Если же применять обычные способы определения этих характеристик практически на "глаз", то это может приводить к погрешностям до 50 %. Предложены методы, которые можно разделить на два класса: лучше идентифицирующие сами глобальные параметры, но хуже - локальную траекторию; и наоборот, хуже идентифицирующие параметры и лучше - траекторию. Естественно, что следует применять те из них, которые более соответствуют конкретной практической задаче. Однако, отметим еще раз, что все эти методы дают как параметры, так и траектории, точнее, чем принятая практика обработки траекторных наблюдений.
Литература
Кацман Ф.М., Музыкантов Г.М., Шмелев А.В. Эксплуатационные испытания морских судов. М.,
Транспорт, 272 с., 1970.
Пашенцев С.В., Мартюк Г.И., Юдин Ю.И. Идентификация параметров математической модели судна.
Программа для ЭВМ. РОСПАТЕНТ, № 2003611468, 19.06.2003. Юдин Ю.И., Пашенцев С.В. Разработка математической модели танкера "Саратов". Деп. ВНИЭРХ, № 1391, РХ-2003.
