ISSN 2075-6836

Ф е д е ра л ь н о е го с уд а р с т в е н н о е б юд ж е т н о е у ч р е ж д е н и е н а у к и

институт космических исследований российской академии наук (ики ран)

Б. Ц. Бахшиян

ОЦенивание и кОррекЦия

параметрОв

движущихся систем

Курс леКций

Механика, управление и инфорМатика МосКва 2012 УДК 519.7 ISSN 2075-6839 ББК 22.2 Б30 Серия «Механика, управление и информатика»

Б. Ц. Бахшиян Оценивание и коррекция параметров движущихся систем Курс лекций Данная книга представляет собой курс лекций по основам теории планирования оптимального эксперимента, который в течение ряда лет читался ведущим научным сотрудником ИКИ РАН Б. Ц. Бахшияном (1944–2011) для аспирантов ИКИ РАН, студентов Московского государственного технического университета им. Баумана и Московского авиационного института. В курсе кратко изложены основные разделы теории оценивания, линейного программирования, планирования эксперимента.

Для преподавателей, студентов и аспирантов факультетов прикладной математики технических вузов, специалистов в области математической теории планирования эксперимента.

Ключевые слова: теория оценивания, гарантирующее оценивание, линейное программирование, симплекс-метод, планирование эксперимента, коррекция движения, линейная импульсная коррекция.

Series “Mechanics, control, and informatics” B. Ts. Bakhshiyan Estimation and correction of parameters of moving systems The book contains lectures on the theory of design of optimal experiment.

This course was delivered by a leading scientist of the Space Research Institute (IKI RAN) B. Ts. Bakhshiyan (1944–2011) for IKI graduates and for students of the Bauman Moscow State Technical University and the Moscow Aviation Institute. Basic concepts of the estimation theory, linear programming and design of experiment are considered.

This book may be of interest to lecturers, graduates and students of applied mathematics departments of technical universities and for specialists in design of experiment.

Keywords: estimation theory, guaranteed estimation, linear programming, simplex method, design of experiment, motion correction, linear impulse correction.

Редактор: Корниленко В. С.

Компьютерная верстка: Комарова Н. Ю.

© Федеральное государственное бюджетное учреждение науки Институт космических исследований Российской академии наук (ИКИ РАН), Предисловие Значительное развитие, которое получили в последние десятилетия методы анализа и планирования эксперимента, диктует необходимость подготовки специалистов, владеющих уже апробированными методами и способных разрабатывать новые.

В связи с этим в ряде вузов читаются курсы, посвященные теории оптимального эксперимента.

Данное пособие написано на основе курса лекций, читавшегося автором на протяжении ряда лет аспирантам Института космических исследований РАН и студентам старших курсов Московского государственного технического университета им. Баумана и Московского авиационного института. Оно ориентировано в первую очередь на студентов факультетов прикладной математики технических вузов и инженеров, активно использующих математические методы.

Пособие не преследует цели дать читателю полное описание современной теории оптимального эксперимента. Оно рассчитано на семестровый курс лекций и содержит основные понятия математической теории оценивания, линейного программирования, постановку и методы решения задач L- и MV-оптимального планирования эксперимента, идеальной линейной импульсной коррекции движения. Излагаются также методы гарантирующего оценивания.

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

Считаем своим долгом выразить благодарность профессору К. В. Семенихину за внимательное прочтение рукописи и ряд существенных замечаний.

Горяинов А. В., Федяев К. С.

Москва, февраль 2012 г.

1. Задача ОЦенивания параметрОв движущейся системы и метОды ее решения 1.1. Модель оценивания Каждый человек в своей жизни и общество в целом ежеминутно решают задачи нахождения оценок (или, короче, оценивания) параметров тех или иных систем. При этом система описывается с помощью совокупности некоторых параметров, характеризующих ее и изменяющихся по каким-либо законам. Бытовые оценки параметров носят обычно описательный характер:

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

Оценивание параметров производится обычно на основе наблюдения за поведением исследуемой системы на некотором интервале времени, при этом измеряются какие-либо величины, связанные с системой. Например, при изучении движения искусственного спутника Земли может измеряться дальность до него в различные моменты времени, и на основе этих наблюдений может быть получена оценка таких параметров движения спутника как скорость или ускорение. Используемая модель в этом случае будет представлять собой систему уравнений, связывающих измеренные значения дальности до спутника с оцениваемыми значениями его скорости и ускорения. Такую модель, связывающую измерения с оцениваемыми параметрами, будем называть далее моделью измерений.

6 1. Задача оценивания параметров движущейся системы и методы ее решения Будем считать, что модель измерений включает m неизвестных параметров. Обозначим через m вектор из этих параметров. Обычно вектор включает начальные условия движения системы и ряд постоянных, определяющих модель движения (например, для движения спутника — параметры атмосферы, аэродинамические характеристики спутника и т. п.).

В соответствии с этим будем полагать, что модель измерений представляется в виде Здесь y — вектор измерений; f : ® — вектор-функция, определяющая модель; — суммарный вектор ошибок модели и измерений.

1.2. несМещенный алгоритМ оценивания Если n = m, то естественно искать оценку вектора как решение системы уравнений В случае, когда число измерений равно числу оцениваемых параметров (n = m) и существует однозначная обратная функция f -1( y), решение последних уравнений имеет вид Если положить = 0, то определенный таким образом оцениватель при n = m обладает свойством =, т. е. оценка параметра совпадает с его истинным значением при отсутствии ошибок исходных данных. Такое свойство будем называть несмещенностью оценивателя.

Рассмотрим теперь понятие несмещенного оценивателя при n > m.

Теорема 1.1 (о несмещенном оценивателе). Пусть — некоторая векторная норма. Оцениватель является несмещенным.

Д о к а з а т е л ь с т в о. Если в (1.1) положить = 0, то минимум в (1.2) достигается при x = — в этом случае норма равна нулю, а при x норма положительна. Поэтому для рассмотренного оценивателя при = 0 имеем =, что и доказывает его несмещенность.

Рассмотрим два вида нормы.

1. Пусть норма евклидова, т. е.

где xi — компоненты вектора x *. Тогда оцениватель (1.2) имеет вид и называется методом наименьших квадратов (МНК). Оценку мнк вектора, полученную с использованием оценивателя (1.3), будем называть оценкой метода наименьших квадратов или, короче, оценкой наименьших квадратов. В более общем случае, когда норма определяется выражением для некоторой положительно определенной матрицы W, которую далее будем называть весовой, оцениватель имеет вид и называется обобщенным методом наименьших квадратов (ОМНК). Соответствующую оценку мнк -W будем называть оценкой обобщенного метода наименьших квадратов или, короче, оценкой ОМНК.

2. Пусть норма имеет вид Тогда оцениватель (1.2) имеет вид и называется методом наименьших модулей (МНМ), а оценка мнм — оценкой наименьших модулей.

* Символом здесь и всюду далее обозначается операция транспонирования.

8 1. Задача оценивания параметров движущейся системы и методы ее решения Качество оценки характеризуется величиной ее ошибки - , а точнее — статистическими характеристиками этой ошибки. Обычно эти характеристики удается найти лишь для линейной модели измерений.

1.3. линеаризация Модели изМерений Пусть 0 — некоторое номинальное значение вектора. Предполагая разность – 0 достаточно малой, проведем линеаризацию модели измерений, оставляя только слагаемые первого порядка. Тогда она запишется в виде где — матрица размерности nm*. Окончательно линеаризованная модель записывается в виде где y = y - f (0 ) + H 0 — приведенный вектор измерений. Эта приближенная модель обычно используется для вычисления оценки = 1, относительно которой модель снова линеаризуется. Процесс последовательных приближений продолжается до тех пор, пока две оценки s и s+1 не окажутся достаточно близs + кими друг к другу. Тогда оценка = принимается за оценку в первоначальной нелинейной модели.

1.4. одноМерная линейная Модель Одномерная линейная модель имеет вид где — неизвестный скалярный параметр.

* То есть матрица, имеющая n строк и m столбцов.

1.5. линейный несмещенный алгоритм и метод наименьших квадратов Можно показать, что оценка наименьших квадратов представляет собой выборочное среднее:

а оценка наименьших модулей — выборочную медиану вектора y.

Из приведенных формул следует, что оценка наименьших квадратов не является устойчивой к аномальным измерениям в отличие от оценки наименьших модулей.

1.5. линейный несМещенный алгоритМ и Метод наиМеньших Квадратов Будем рассматривать только линейную модель оценивания. Нас будет интересовать значение некоторого скалярного параметра, являющегося линейной функцией вектора :

(далее параметр l будем называть контролируемым), а также совокупности подобных параметров где L s — вектор контролируемых параметров; B — матрица размерности ms, b и B известны. Если по результатам измерений найдена оценка вектора, то оценки указанных параметров можно найти по формулам Для линейной модели оценка ОМНК, соответствующая весовой матрице W, имеет вид где Как видим, если выполняется условие det( HWH ) 0, то ОМНК представляет собой линейный относительно вектора y несмещенный алгоритм (доказательство несмещенности предоставляется читателю). При этом 10 1. Задача оценивания параметров движущейся системы и методы ее решения Отвлечемся от метода наименьших квадратов и рассмотрим сразу линейную оценку вектора L:

где X — матрица размерности ns, определяющая эту оценку.

Условие несмещенности линейного алгоритма имеет вид Тогда ошибка оценки будет равна Таким образом, ОМНК дает линейную несмещенную оценку любого векторного параметра L = B (а значит, и любого скалярного параметра l = b ). Точнее, каждой весовой матрице W однозначно соответствует одна матрица X для линейного оценивания вектора L = B (или один вектор x для оценивания параметра l = b ).

Рассмотрим соотношение между оценкой ОМНК и линейной несмещенной оценкой с другой точки зрения. Пусть известна матрица X размерности ns, удовлетворяющая условию несмещенности (в скалярном виде это условие имеет упомянутый выше вид Hx = b). Как найти весовую матрицу W такую, чтобы оценка ОМНК вектора L = C, соответствующая этой матрице, имела вид L = X y ? Эта задача представляет собой задачу нахождения неотрицательно определенной матрицы W, являющейся решением матричного уравнения Такая задача рассматривалась различными авторами (см., например, [Rao, 1975]) и, более подробно, нами [Бахшиян, 1983].

Существует целое множество весовых матриц, для каждой из которых оценка наименьших квадратов совпадает с заданной несмещенной оценкой. Явное описание этого множества здесь не приводится.

Таким образом, на основе предыдущих рассуждений можно сделать вывод о том, что множество всех линейных несмещенвычисление точности оценок при известной ковариационной матрице ошибок… ных оценок любого векторного параметра L = B совпадает с множеством всех оценок ОМНК, получающихся при всевозможных весовых матрицах. Поэтому в дальнейшем будем анализировать только линейные несмещенные алгоритмы.

Замечание 1.1. Рассмотренный выше метод наименьших модулей не является линейным несмещенным оценивателем. Поэтому он, вообще говоря, не эквивалентен никакой оценке ОМНК.

1.6. вычисление точности оценоК При известной Ковариационной Матрице ошибоК изМерений.

теореМа гаусса – МарКова Пусть вектор ошибок представляет собой случайный вектор с математическим ожиданием E(), равным нулевому вектору, и известной матрицей ковариаций K E ( ). Тогда дисперсия ошибки оценки l = x y скалярного параметра l = b равна где l = l - l. Данная формула является универсальной, так как она справедлива для любого оценивателя x и используется для вычисления не только дисперсии оценок скалярных параметров, но и для определения ковариационной матрицы оценок векторных параметров. В этом случае возникает вопрос о нахождении оптимального оценивателя, для которого дисперсия оценки будет минимальной. Можно показать, что минимум дисперсии по оценивателю x при условии несмещенности Hx = b достигается при Это устанавливается с помощью метода множителей Лагранжа. Можно заметить, что оцениватель получился такой же, как если бы использовался ОМНК с весовой матрицей W = K –1. Таким образом, с точки зрения минимизации дисперсии оценки любого скалярного параметра оптимальным является ОМНК с весовой матрицей, обратной матрице ковариаций ошибок измерений. Этот известный результат называется теоремой Гаусса – Маркова. Оптимальная дисперсия при этом равна 12 1. Задача оценивания параметров движущейся системы и методы ее решения где — ковариационная матрица оценки вектора.

Замечание 1.2. Можно показать, что матрица является неотрицательно определенной, если K — ковариационная матрица оценки, полученной любым линейным несмещенным алгоритмом.

1.7. ЭКвивалентность Множеств всевозМожных линейных несМещенных оценоК При различноМ выборе Мешающих ПараМетров Рассмотрим снова линейную модель измерений где y — вектор измерений с компонентами yi, i = 1, …, n;

— матрица размерности mn, которую будем называть матриm цей плана измерений; — вектор оцениваемых параметров;

— суммарный вектор ошибок измерений. Пусть — скалярный контролируемый параметр; — ошибка модели (1.6). Такая модель на практике получается линеаризацией некоторой реальной нелинейной модели, при этом представляет собой отклонение вектора некоторых неизвестных параметров относительно их априорных значений. Поэтому можно считать, что априорное значение вектора также известно и равно нулю.

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

В соответствии с изложенным, к уравнениям модели (1.4) следует добавить еще m уравнений 1.7. Эквивалентность множеств всевозможных линейных несмещенных оценок… На практике обычно некоторый подвектор II вектора принимают равным своему априорному значению II (II ), а оцениваемым вектором считают совокупность I остальных компонент вектора. Это обусловлено, в частности, следующими причинами:

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

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

При указанном подходе компоненты вектора II называют мешающими параметрами. В связи с разбиением вектора запишем известные вектор и матрицу в (1.4), (1.6) в составном виде:

Для модели с мешающими параметрами обобщенный метод наименьших квадратов записывается следующим образом:

где — вектор невязок, рассчитанный при априорно заданном векторе II, W — весовая матрица. При этом оценка интересующего нас параметра имеет вид Минимизация в (1.9) приводит к известному результату, согласно которому оценка метода наименьших квадратов есть линейная функция вектора измерений и априорного значения вектора мешающих параметров. А именно, решение оптимизационной задачи (1.9), согласно п. 1.5, имеет вид 14 1. Задача оценивания параметров движущейся системы и методы ее решения Подставляя (1.12) в (1.11), получаем выражение для оценки (1.11) в зависимости от исходных данных:

где — вектор коэффициентов, называемый в дальнейшем линейным оценивателем.

Таким образом, оценка методом наименьших квадратов (1.9) приводит для линейной модели к линейной оценке относительно исходных данных y и II. Прямыми вычислениями найдем ошибку l оценки (1.13). Для этого подставим в (1.13) выражение (1.4) для вектора измерений и вычтем из полученного выражения истинное значение (1.6) оцениваемого параметра.

Имеем

II II II II I I II

Согласно (1.14), Поэтому окончательно получаем где — ошибка априорного значения вектора мешающих параметров или, как будем далее говорить, ошибка априорного измерения.

Согласно (1.17), ошибка оценки равна нулю, если равны нулю ошибки исходных данных, II. Поэтому алгоритм метода наименьших квадратов является несмещенным. Несмещенным будет также любой линейный алгоритм (1.13) с произвольным оценивателем, удовлетворяющим условию (1.14). В связи с этим условие (1.14) называется условием несмещенности.

Таким образом, ошибка линейной оценки (1.13) для модели (1.6) имеет вид (1.17).

В заключение параграфа приведем еще одну форму записи выражения (1.17). Введем расширенный вектор измерений 1.7. Эквивалентность множеств всевозможных линейных несмещенных оценок… включающий наряду с измерениями и априорные значения вектора мешающих параметров, а также расширенный вектор ошибок соответствующий вектору (1.19). Напомним, что в вектор измерений y могут входить априорные значения компонент вектора определяемых параметров II. Введем также расширенную матрицу плана измерений, соответствующую вектору (1.19) где символами 0 и I обозначены матрица с нулевыми элементами и единичная матрица соответствующих размерностей*.

При введенных обозначениях и условии, что вектор мешающих параметров полагается равным своему априорному значению, модель измерений (1.4) записывается в виде а линейный алгоритм оценивания (1.13) — в виде где составной вектор представляет собой обобщенный оцениватель. Этот оцениватель удовлетворяет условию несмещенности которое следует из аналогичного условия (1.16). Подставляя (1.22) в (1.23) и вычитая из полученного выражения истинное значение l = b, получаем с учетом условия (1.25) выражение для ошибки оценивания:

* Далее также всюду будем обозначать символами 0n и 0mn соответственно вектор размерности n и матрицу размерности mn, состоящие из нулевых элементов; символом In — единичную матрицу размерности n.

16 1. Задача оценивания параметров движущейся системы и методы ее решения Отметим, что это выражение есть просто другая запись равенства (1.17) при введенных обозначениях.

Итак, оценка наименьших квадратов с любой весовой матрицей является линейной несмещенной оценкой общего вида (1.23), (1.25). Может быть доказано и обратное утверждение:

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

Таким образом, имеет место совпадение множества линейных несмещенных оценок (1.23) (т. е. всех оценок, определяемых из (1.23) при всевозможных x, удовлетворяющих условию (1.25)), и множества оценок методом наименьших квадратов, соответствующих всевозможным весовым матрицам при произвольном выборе вектора II мешающих параметров (II II ).

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

2. ОснОвные сведения иЗ теОрии линейнОгО прОграммирОвания 2.1. ПостановКа задачи и сиМПлеКс-Метод Будем рассматривать задачу линейного программирования (ЛП) следующего вида:

где x1, …, xn — аргументы линейной функции z(x), которую часто называют целевой; c1, …, cn — числовые коэффициенты, a1, …, an; b.

В векторно-матричной форме эта задача имеет вид Здесь x n ; c n ; A — матрица коэффициентов ограничений размерности mn; b m — вектор правых частей ограничений.

В дальнейшем будем полагать n m и rang A = m (невыполнение условия n m может привести к неразрешимости системы ограничений, а условие rang A < m говорит о наличии избыточных ограничений либо об их несовместности).

Выделим в матрице A невырожденную подматрицу B размерности mm, которую будем называть базисной матрицей или, короче, базисом. Перенумеровав при необходимости переменные, приведем матрицу A к виду где N — матрица размерности m(n – m). Соответственно разобьем векторы x и c:

18 2. основные сведения из теории линейного программирования Вектор x = ( x B, 0,, 0), такой что Ax = b и x 0, будем называть допустимым базисным вектором задачи ЛП, соответствующим базису B, который, в свою очередь, будем называть допустимым. Вектор x, такой что c x = z, будем называть оптимальным вектором или решением задачи (2.1). Можно доказать, что среди оптимальных векторов задачи (2.1) всегда есть базисный. Базис, соответствующий оптимальному вектору, будем также называть оптимальным.

Таким образом, поиск решения задачи (2.1) может быть сведен к перебору всевозможных базисных векторов и выбору среди них того, для которого значение целевой функции z(x) минимально. Идея симплекс-метода состоит в том, чтобы сделать такой перебор направленным, рассматривая на каждом последующем шаге только те базисные векторы, на которых значение целевой функции не больше, чем на текущем.

Представим вектор xB в виде Полагая xN = 0, получим выражение для вектора xB, соответствующего базисному вектору:

Введем вспомогательный вектор из условия Тогда справедливо равенство cB x B = cB B b = b.

Пусть x — текущий допустимый базисный вектор, соответствующий базису B. Рассмотрим произвольный вектор x x и найдем z( x ) :

где = c - A. Отсюда следует, что текущий базис B оптимален тогда и только тогда, когда В частности, достаточным условием оптимальности является выполнение неравенства Предположим, что это достаточное условие не выполняется, т. е. существует номер s такой, что s < 0. Увеличивая от нуля значение компоненты xs вектора x, как следует из формулы (2.4), будем уменьшать значение целевой функции. При этом, однако, должно выполняться условие откуда или где = B as. Из последнего неравенства находим — максимально возможное значение x s, при котором вектор x остается допустимым. Если же i 0 при всех i = 1, …, m то компоненту xs можно неограниченно увеличивать, не нарушая при этом условий допустимости. Из формулы (30) следует, что в этом случае целевая функция z( x ) не ограничена. В остальных случаях, полагая xs =, получим x На основе проведенных рассуждений можно построить следующий итерационный алгоритм нахождения решения x задачи ЛП.

1. Пусть B — некоторый базис, допустимый для задачи ЛП.

Вычисляется вспомогательный вектор из условия 2. Вычисляется вектор симплекс-разностей :

Заметим, что условие (2.6) эквивалентно условию поэтому на практике при заданном базисе B вычисляют только вектор N = ( m+1,..., n ).

20 2. основные сведения из теории линейного программирования 3. Вычисляется минимальная симплекс-разность 4. Проверяется достаточное условие оптимальности Если это условие выполняется, то текущий базис B оптимален, и вычисления завершаются. В противном случае столбец as вводится в базис.

5. Вычисляется вектор координат вектора as в базисе B, и ищется величина Если вектор не содержит положительных компонент, то целевая функция z(x) не ограничена на множестве допустимых векторов, и вычисления завершаются. В противном случае столбец ar выводится из базиса.

6. Производится пересчет значений базисных переменных, целевой функции и матрицы B по формулам где gij — элементы матрицы B ; gij — элементы матрицы B.

После этого происходит возврат к шагу 1. Новый допустимый базис B будет состоять из векторов ai, i = 1, …, m, i r и вектора as.

Легко видно, что если на каждой итерации симплекс-метода величина положительна, то приведенный алгоритм за конечное число итераций позволяет либо найти решение задачи, либо установить, что целевая функция не ограничена.

2.2. оценКа близости теКущего решения К оПтиМальноМу для неКоторых тиПов задач 1. Случай, когда ci = 1 для всех i. Этот случай имеет место, например, для задачи ЛП, к которой сводится C-задача планирования эксперимента, которая будет рассмотрена позднее. Из формулы (2.1) получаем где x — решение задачи (2.1). Отсюда получаем оценку для оптимума:

2. Случай наличия ограничения вида В этом случае из (2.4) легко получить оценку 2.3. Построение доПустиМого базисного решения В приведенных выше рассуждениях подразумевалось, что в матрице A всегда можно указать некоторый начальный допустимый базис B. Однако при решении практических задач непосредственно сделать это бывает невозможно. Для построения допустимого базиса задачи ЛП применяются специальные методы, один из которых описан ниже. Для простоты изложения будем далее полагать, что b 0, чего легко добиться изменением знаков строк матрицы условий.

Рассмотрим следующую вспомогательную задачу:

Здесь e — вектор с единичными компонентами; y — вектор вспомогательных переменных, называемых искусственными; I — единичная матрица размерности mm.

22 2. основные сведения из теории линейного программирования Задача (2.15) представляет собой задачу линейного программирования с матрицей условий размерности m(n + m). Однако для нее можно непосредственно указать допустимый базис, состоящий из столбцов матрицы I.

Если значение целевой функции w на оптимальном векторе задачи (2.15) положительно, то построение допустимого базиса исходной задачи ЛП невозможно, так как в задаче (2.15) такому базису соответствовало бы нулевое значение целевой функции w.

Если же значение целевой функции w на оптимальном векторе равно нулю, и оптимальный базис не содержит столбцов матрицы I, то этот базис, очевидно, является допустимым для исходной задачи ЛП. Если же w = 0 и оптимальный базис задачи (2.15) содержит столбцы матрицы I, то соответствующие базисные компоненты вектора y будут нулевыми. Эти столбцы необходимо заменить на столбцы матрицы A. Построенный таким образом вырожденный базис будет допустимым для исходной задачи ЛП.

2.4. геоМетричесКий сПособ решения двуМерной задачи линейного ПрограММирования сПециального вида Рассмотрим двумерную задачу ЛП вида в которой целевая функция имеет единичные коэффициенты, a1, …, an, b.

Решение такой задачи может быть построено геометрически (рис. 1). Для этого на координатной плоскости из начала координат строятся векторы a1, …, an, b. Затем выбирается базис из каких-либо двух векторов ai (на рисунке начальный базис состоит из векторов a1 и a3). Потом проводится прямая через концы векторов выбранного базиса. Если среди векторов a1, …, an нет векторов, пересекающих проведенную прямую, то базис оптимален. На рис. 1 вектор a2 пересекает эту прямую.

Если его ввести в базис вместо вектора a3, то, как нетрудно заметить, коэффициенты разложения вектора b по векторам a1, a 2.4. геометрический способ решения двумерной задачи линейного программирования… также будут положительными, и сумма этих коэффициентов будет меньше, чем сумма коэффициентов на предыдущем шаге, т. е. значение целевой функции уменьшится. Далее снова проводится прямая через концы векторов a1 и a2 и т. д. Такой упорядоченный перебор базисов и составляет идею симплексного метода.

3. классические Задачи планирОвания эксперимента и их сведение к Задачам ОптимальнОй идеальнОй импульснОй кОррекЦии параметрОв системы 3.1. усреднение Модели изМерений Пусть m — вектор неизвестных параметров системы, и при каждом значении i = 1, …, n проводится ri измерений yij, j = 1, …, ri некоторой функции H i. Обычно считают, что измерения yij при заданном i проводятся в заданный момент времени ti (возможен и другой параметр привязки всех измерений yij при заданном i, например дальность полета при движении точки по траектории). Тогда можно измерения в этот момент считать сеансом измерения. Например, космические измерения проводятся сеансами, соответствующими небольшим интервалам времени, в течение которых есть радиовидимость космического объекта. Эти интервалы часто можно при планировании космического эксперимента заменять отдельными моментами ti.

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

Согласно теореме Гаусса – Маркова, в классе линейных несмещенных алгоритмов оценивания минимальную дисперсию оценки любого параметра l = b дает МНК, соответствующий весовой матрице W = I. Учитывая это, осредним модель оценивания.

Лемма 3.1. При использовании МНК исходная модель измерений эквивалентна усредненной модели где — средние арифметические указанных выше ri измеренных значений yij и ошибок ij. Эквивалентность данных моделей означает, что оценки параметра, получаемые при использовании этих моделей, совпадают. При этом усредненные ошибки i являются некоррелированными между собой случайными величинами, математические ожидания которых равны нулю, а дисперсии равны Для доказательства леммы согласно теореме Гаусса – Маркова запишем дисперсию оценки МНК некоторого скалярного параметра l = b, соответствующей весовой матрице W = I, как минимальную дисперсию на классе всех линейных несмещенных алгоритмов. Получаем где коэффициенты xij оценивателя находятся из условия минимума дисперсии оценки при условиях несмещенности:

Минимум по xij суммы в (3.2) при условии (3.3) и фиксиj =1 x рованных xi достигается при xij = i для всех j, т. е. все повторяri ющиеся измерения в оптимальном случае оцениваются с одним оценивателем (это ясно и из соображений симметрии). Тогда (3.2) переписывается в виде 26 3. классические задачи планирования эксперимента… Согласно теореме Гаусса – Маркова равенство (3.4) доказывает лемму 3.1, так как линейный оцениватель с коэффициентами xi соответствует модели измерений (3.1).

Замечание 3.1. Так как в (3.4) суммирование записано по всем i, то, чтобы исключить суммирование по тем элементам, где ri = 0, целевая функция должна быть доопределена следующим образом:

3.2. ПостановКа задачи Пусть задано общее число измерений:

Вместо чисел ri будем использовать вектор p с координатами Далее будем, как это обычно делают, пренебрегать целочисленностью ri и считать, что p — непрерывный план эксперимента, т. е.

любой вектор из симплекса контролируемые параметры (обычно s m ). При заданном плане p рассмотрим линейные оценки каждого из этих параметров:

где числа ij, i = 1, …, n определяют линейный оцениватель параметра lj. Этот оцениватель дает несмещенные оценки парамеПостановка задачи тров b тогда и только тогда, когда он удовлетворяет условиям несмещенности т. е. bj есть линейная комбинация тех векторов Hi, для которых pi > 0. Далее будем обозначать символом n множество тех p из n, для которых уравнения (3.6) разрешимы относительно ij.

Для каждого p n рассмотрим наилучшие линейные несмещенные оценки, для которых оцениватель определяется из условия минимальных дисперсий при линейном несмещенном оценивании:

Задача оптимального планирования эксперимента состоит в нахождении плана p n, доставляющего нижнюю грань заданному критерию оптимальности L(p) [Бахшиян, 1970, 1983, 1989; Бахшиян, Соловьев, 1998; Бахшиян и др., 2000]:

Будем рассматривать два критерия оптимальности, предполагая далее, что оценки (3.5) являются наилучшими нелинейными несмещенными. Первый критерий называется критерием L-оптимальности. Второй критерий будем называть MVs-критерием (при s = m и lj = j, j = 1, …, m критерий (3.10) называют просто MV-критерием). В зависимости от критерия (3.9) или (3.10) задачу (3.8) будем соответственно называть L-задачей или MVs-задачей.

При решении задачи (3.8) будем предполагать, что 28 3. классические задачи планирования эксперимента… (т. е. ранг составной матрицы максимален), что равносильно допущению о возможности линейного несмещенного оценивания всех компонент вектора в случае использования всех n групп измерений.

Замечание 3.2. Нетрудно показать, что суммирование по i p в соотношениях (3.5)–(3.7) можно заменить суммированием по i = 1, …, n, если доопределить слагаемые в (3.7) при pi = следующим образом:

В дальнейшем будем считать указанные замены при суммировании и оптимизации произведенными.

Замечание 3.3. Согласно теореме Гаусса – Маркова при выполнении условия (3.6) наилучшая нелинейная несмещенная оценка совпадает с оценкой взвешенного метода наименьших квадратов l = b, где — любое решение нормальных уравнений — нормированная на число N информационная матрица. Эти уравнения всегда совместны (доказательство этого факта предоставляется читателю), но решение будет единственно только в случае невырожденной матрицы M(p). План p в этом случае называется невырожденным. В противном случае план p является вырожденным, и решение уравнений (3.12) неоднозначно. Однако условие (3.6) эквивалентно тому, что оценка b каждой функции b, j = 1, …, s одна и та же для всех решений уравнения (3.12).

При этом параметры b называются оцениваемыми. Таким обj разом, условия несмещенности и оцениваемости эквивалентны.

Далее, условия несмещенности (3.6) можно трактовать как принадлежность векторов bj линейному многообразию, порожденному векторами p1 H1,, pn H n или, эквивалентно, столбцами информационной матрицы. Это означает, что найдутся векторы 3.3. скалярная задача планирования эксперимента и алгоритм ее решения j m, такие что выполняется условие оцениваемости параметров b :

3.3. сКалярная задача Планирования ЭКсПериМента и алгоритМ ее решения При s = 1 L-задача (а также MVs-задача) сводится к минимизации дисперсии Dl, где l = l1 и называется C-задачей планирования эксперимента. Согласно (3.4) C-задача записывается в виде Проведем далее минимизацию дисперсии по pi, пренебрегая тем, что ri должно быть целым. Это оправдано при достаточно большом N. Получим, используя метод множителей Лагранжа, что оптимальный план измерений есть При этом для заданного оценивателя с коэффициентами xi Замечание 3.4. Указанная оптимизация по ri здесь проведена не совсем строго (так же как и в основополагающей работе [Elfving, 1952]). Действительно, мы проводили оптимизацию по pi только при условии p1 + … + pn = 1. Однако условие несмещенности может выполняться не при всех таких pi, т. е. оценка параметра l может быть получена не всегда. В связи с этим условие несмещенности также называют условием оцениваемости. Поэтому множество допустимых pi может быть несколько меньше, а искомый минимум — не меньше. Однако ввиду того, что получившееся решение приводит к допустимым оценкам, оно является верным.

30 3. классические задачи планирования эксперимента… Оптимизация величины D0, полученной при оптимальных ri, по оценивателю x приводит к задаче линейного программирования Введем новые переменные zi = xi, zi 0. Тогда xi Hi = zi ai, где При этом оптимальную задачу можно записать в виде где i = {H i, - H i } — множества, состоящие из двух векторов (рис. 2).

Сформулированную задачу можно интерпретировать следующим образом: минимизировать сумму коэффициентов zi > разложения вектора b по векторам ±Hi. Эта задача представляет собой задачу линейного программирования, геометрический Рис. 2. Построение выпуклой оболочки множеств i 3.4. оценивание параметров параболической траектории по измерениям дальности способ ее решения (в двумерном случае) обсуждался в предыдущей главе. У этой задачи существует решение, содержащее не более m отличных от нуля коэффициентов, т. е. вектор b лежит внутри многогранного угла, образованного некоторыми m векторами ±Hi (базисом), причем гиперплоскость, проведенная через их концы, отсекает от вектора b наименьшую (считая от конца вектора) часть по сравнению с другими комбинациями из m векторов. Легко заметить, что указанная гиперплоскость, проходящая через концы векторов оптимального базиса, принадлежит выпуклой оболочке всех 2n векторов ±Hi, i = 1, …, n (см. рис. 2).

Таким образом, можно сделать следующий важный вывод.

При решении скалярной задачи планирования эксперимента следует выбрать оптимальные m измерений и найти оцениватель для усредненной модели из соотношения где H — матрица размерности mm, составленная из векторов Hi, входящих в оптимальный базис. Это соответствует тому, что вектор находится из m уравнений где нумерация по i условна. После этого оптимальные значения величин ri находятся из условия их пропорциональности величинам xi.

3.4. оценивание ПараМетров ПараболичесКой траеКтории По изМеренияМ дальности Рассмотрим в качестве примера плоское движение материальной точки, брошенной под углом к горизонту с начальной скоростью v0. Введем следующие обозначения: r, v — векторы положения и скорости точки в момент времени t, r0, v0 — векторы положения и скорости в начальный момент времени t0 = 0, Движение точки описывается моделью 32 3. классические задачи планирования эксперимента… где ма на максимальную высоту. Пусть на интервале [0, tk] могут производиться измерения дальности l = r (t ) до движущейся точки. Ошибки дальности будем считать случайными величинами с нулевым математическим ожиданием и стандартными отклонениями l = 1 м.

В качестве вектора оцениваемых параметров рассмотрим начальный вектор состояния = X(0). При этом будем считать, что положение точки в начальный момент известно, а именно x(0) = 0, y(0) = 0, т. е. число оцениваемых параметров равно двум (оцениваются компоненты вектора v0). В качестве контролируемого параметра выберем дальность до точки в момент падения.

Требуется:

• решить задачу оптимального оценивания дальности на интервале [0, tk], т. е. найти оптимальные значения моментов измерений t1, t2 и соответствующую среднеквадратичную • доказать, что при любом tk справедливо t2 = tk;

• интерпретировать результат как C-задачу планирования эксперимента, найти соответствующие доли p1, p2 общего числа измерений в моменты времени t1 и t2.

Запишем линеаризованную модель измерений в виде где H ( H r, H v ) — составной вектор-строка.

Нахождение оптимальных моментов измерений, согласно изложенной выше теории, сводится в данном случае к решению двумерной задачи линейного программирования (используется только вектор Hv(t)) 3.4. оценивание параметров параболической траектории по измерениям дальности где b = (T, 0) (переменные, по которым производится оптимизация, обозначены здесь через z, так как обычное обозначение x уже используется для первой компоненты вектора r (t)). Будем искать решение этой задачи геометрически в соответствии с подходом, изложенным в предыдущем разделе. При этом для простоты и наглядности будем пренебрегать дискретностью моментов измерений, рассматривая вместо моментов ti непрерывную переменную t.

Компоненты матрицы Hv(t) вычисляются из соотношений Производные от конкретных измеряемых функций вычисляются по формулам Окончательно векторы Hv(t) для каждого t будут иметь вид Годограф функции Hv(t) имеет вид, показанный на рис. 3.

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

Построение линейной оболочки и оптимальный базис показаны на рис. 4. Нетрудно заметить, что при любом значении tk [0, T ] вектор Hv(tk) входит в состав оптимального базиса, т. е. выполняется условие t2 = tk.

Момент времени t1 может быть найден из условия касания годографа и прямой, проходящей через его крайнюю точку, соответствующую моменту времени tk. После этого из решения 34 3. классические задачи планирования эксперимента… задачи (3.17) определяются z1 и z2, а по ним — p1 и p2. Значение среднеквадратической ошибки, согласно (3.16), равно оптимальному значению целевой функции задачи (3.17).

Возьмем, например, tk = 80 с. Для этого случая найдем:

t2 = tk = 80 с, t1 = 38 с, x1 = –1,49, x2 = 1,96, p1 = 0,43, p2 = 0,57, 1 = 3,46.

Рис. 4. Построение оптимального базиса в задаче (3.17) 3.5. сведение L-задачи к задаче оптимальной линейной импульсной коррекции 3.5. сведение L-задачи К задаче оПтиМальной линейной иМПульсной КорреКции Рассмотрим теперь L-задачу с критерием оптимальности (3.9). В соответствии с (3.7) целевую функцию (3.9) можно преобразовать следующим образом:

Здесь вектор ui новых переменных и его норма определяются из условия выражение понимается так же, как и в замечании 3.2, и использованы обозначения:

т. е. b, Ui — матрица размера Ms, где M = ms, 0m — нулевой вектор размерности m.

Теорема 3.1. Оптимальные значения вектора p, линейного оценивателя u, а также значение L-задачи находятся из соотношений где ui — решение задачи 36 3. классические задачи планирования эксперимента… Д о к а з а т е л ь с т в о. Подставляя по-новому записанную целевую функцию в (3.8) и меняя затем порядок оптимизации по p и ui, найдем аналитически экстремум по вектору p при помощи метода множителей Лагранжа так же, как при решении C-задачи планирования эксперимента. При фиксированных переменных ui этот минимум достигается при откуда и получаем указанный результат.

Замечание 3.5. При s = 1 задача (3.20) является задачей линейного программирования. При этом симплекс-методом находится базисное решение, для которого отличны от нуля не более m переменных pi.

Задача (3.20) представляет собой задачу оптимальной линейной импульсной коррекции, где ui и Ui — соответственно корректирующий импульс и матрица его влияния в i-й коррекции, вектор b — требуемая суммарная величина коррекции, а целевая функция характеризует суммарные затраты на коррекцию. Подробнее эта задача и алгоритм ее решения рассматриваются в следующем разделе.

3.6. сведение MVs-задачи К ПараМетричесКой задаче оПтиМальной линейной иМПульсной КорреКции и возМожный алгоритМ ее решения 3.6.1. Сведение MVs-задачи к проектной задаче коррекции Критерий (3.10) можно записать в виде где искусственно вводимый вектор переменных µ = (µ1,,µ s ) принадлежит симплексу Таким образом, задача (3.8) становится минимаксной.

Можно доказать следующий результат.

Теорема 3.2. Если параметры b, j = 1, …, s являются оцениj ваемыми в смысле замечания 3.3, то в минимаксной задаче (3.21) можно поменять местами операции минимума и максимума.

Используя утверждение теоремы 3.2 и записывая вариационное представление дисперсий (3.7), получим:

Условие j > 0 вводится здесь, чтобы исключить неопределенности вида [0 ]. Переставляя операции минимизации по pi и ij, меняя порядок суммирования по i и j и затем проводя оптимизацию по pi, как и при доказательстве теоремы 3.1, находим где L() есть значение подзадачи Если при этом µ j — оптимальное значение в (3.22), а ij — оптимальные значения в (3.23), то Будем рассматривать пока только векторы с положительными компонентами (общий случай будет рассмотрен ниже). При > 0 можно по аналогии с (3.18) ввести новые переменные Тогда подзадача (3.23) запишется в виде 38 3. классические задачи планирования эксперимента… Выражение (3.22) с учетом (3.24) представляет собой проектную задачу оптимальной линейной импульсной коррекции*, в которой оптимальный корректируемый вектор b(µ ) отыскивается из условия, что минимальные затраты на коррекцию при каждом µ будут максимальными среди всех µ s.

Замечание 3.6. Подзадача (3.24) может быть записана в виде оптимальной задачи коррекции (3.20) с заменой L на L() и b на b(). Поэтому она при каждом фиксированном эффективно решается методом генерации столбцов, который будет рассмотрен далее.

3.6.2. редукция задачи (3.22) к задаче МногоМерной МакСиМизации и алгоритМ ее решения Нетрудно показать, что задача, двойственная к подзадаче (3.24), имеет вид:

где Замечание 3.7. Задача (3.25) представляет собой обобщенную задачу линейного программирования, значение которой, как известно [Данциг, 1966, Лэсдон, 1975], непрерывно зависит от вектора b() и, следовательно, от вектора µ 0. Поэтому соотношение двойственности (3.25) имеет место для всех µ s. Если при этом некоторые компоненты вектора равны нулю, то соответствующие уравнения в условиях задачи (3.24) не учитываются.

Необходимо только следить, чтобы соответствующие векторы bj удовлетворяли условию (3.6).

* О проектной задаче оптимальной линейной импульсной коррекции см. ниже, разд. 4.2, формула (4.5).

Соотношение (3.25) позволяет записать (3.22) в виде задачи максимизации по и :

Используя симметричность эллипсоидов i, нетрудно показать, что максимум по при фиксированном достигается при где j определяются из условия = (1,, ). Этот максимум равен Теперь в соответствии с (3.27) можно предложить конструктивный способ покоординатного увеличения целевой функции в (3.26). Возьмем некоторое значение µ, µ > 0 и опредеk k лим последовательности, следующим образом. Вектор k + 1 находится по k из формул (3.27). Вектор k определяется по k из формулы kB = (1, …, 1), где оптимальная базисная матрица B находится, согласно замечанию 3.6, путем решения задачи (3.24) методом генерации столбцов. Согласно (3.26), получаем Замечание 3.8. Процедура (3.28) позволяет монотонно увеличивать целевую функцию в (3.26). Однако, если некоторые компоненты вектора обращаются в нуль, то, как показывает опыт, это может не привести к нахождению глобального максимума в (3.26). В работе [Войсковский, 1998] MV-задача сводится к обобщенной задаче линейного программирования, что позволяет эффективно находить глобальный оптимум. В то же время можно показать, что, если в процедуре (3.28) существует предел µ > последовательности, то он равен оптимальному значению µ в задаче (3.25).

4. Задача ОптимальнОй линейнОй импульснОй кОррекЦии и алгОритм ее решения 4.1. задача оПтиМальной КорреКции При известноМ КорреКтируеМоМ веКторе Задача оптимальной коррекции при известном векторе корректируемых параметров имеет большой самостоятельный интерес и разнообразные практические приложения. В данной работе эта задача изучается ввиду того, что к ней сводится L-задача планирования эксперимента. Пусть l s — некоторый вектор параметров системы, который может быть изменен путем импульсной коррекции ее траектории. Такая коррекция заключается в мгновенном изменении некоторых фазовых координат в моменты времени t1, …, tn. Изменение в каждый момент t1 будем характеризовать корректирующим вектором (импульсом) ui разs мерности si, принадлежащим евклидовому пространству i.

Например, в случае коррекции траектории космического аппарата, если импульс ui может производиться в произвольном направлении, то si = 3, i = 3 — трехмерное евклидово проs странство; если же заранее задана плоскость или прямая, которым должен принадлежать импульс, то si равно двум или единице, а i есть плоскость или прямая. Будем рассматривать идеальную линейную коррекцию, т. е. предполагать, что импульсы производятся без ошибок и изменение вектора l в момент ti под действием импульса ui равно Ui ui, где Ui — матрица размерности ssi, характеризующая влияние импульса и называемая матрицей влияния. Тогда для изменения вектора l на векторную величину b s нужно импульсы ui произвольно выбирать из условия 4.1. Задача оптимальной коррекции при известном корректируемом векторе Пусть затраты на коррекцию в момент ti характеризуются величиной pi (ui), где pi (·) — некоторая норма в i. Например, в случае коррекции траектории космического аппарата путем ориентации в пространстве с использованием одного двигателя, жестко связанного с аппаратом, затраты на коррекцию малого импульса пропорциональны величине Если же для коррекции используются двигатели, действующие вдоль трех осей координат, то где uij — составляющие вектора ui вдоль направления действия двигателей.

При сделанных допущениях поставим задачу оптимизации суммарных затрат на коррекцию. Эта задача имеет вид Решение данной задачи определяет и оптимальные моменты приложения импульса, так как, если ui = 0, то в момент ti импульс не проводится. Представим каждый вектор ui в виде где вектор i i принадлежит множеству { i : pi ( i ) =1}. Тем самым, xi есть длина вектора ui в метрике, определяемой нормой pi(·), а i есть единичный вектор в этой метрике, направленный вдоль импульса.

Используя введенные новые переменные, запишем оптимальную задачу коррекции в виде 42 4. Задача оптимальной линейной импульсной коррекции и алгоритм ее решения где i = {ai : ai = U i i, pi ( i ) =1, i i }. Эта задача по виду напоминает задачу линейного программирования, однако здесь переменными оптимизации являются не только числа xi, но и векторы ai в линейных условиях-равенствах, выбираемые из множеств i. Такая задача называется обобщенной задачей линейного программирования [Данциг, 1966, Бахшиян и др., 1980].

Ее можно также назвать многопараметрической задачей линейного программирования. Приведение задачи коррекции к обобщенной задаче линейного программирования и алгоритм ее решения были предложены М. Л. Лидовым и описаны в работе [Лидов, 1971]. Фактически этот алгоритм представляет собой реализацию так называемого метода генерации столбцов в обобщенном линейном программировании. В этой же статье, а также в книге [Бахшиян и др., 1980] рассмотрены более общие постановки задач коррекции. Далее изложение следует этим работам.

Геометрический способ решения обобщенной задачи линейного программирования аналогичен рассмотренному выше геометрическому способу решения обычной задачи линейного программирования. Для случая, когда нормы pi являются евклидовыми, множества i представляют собой линейное преобразование сфер i =1, т. е. являются эллипсоидами. Такой случай изображен на рис. 5. Если вектор b находится внутри многогранного угла, образованного некоторыми k векторами ai i, i = 1, …, k, то векторы ai, i = 1, …, k дают решение задачи в следующем смысле (см. рис. 5).

4.1. Задача оптимальной коррекции при известном корректируемом векторе Нужно провести k s корректирующих импульсов, величина каждого импульса равна xi, направление определяется вектором i. При этом случай k < s уже не является исключительным. Например, если вектор b совпадает по направлению с вектором ai, принадлежащим выпуклой оболочке i, то k = 1.

Симплексный алгоритм также сводится к вводу в базис вектора as, который пересекает гиперплоскость, проходящую через концы векторов базиса. Проверка достаточного условия оптимальности, таким образом, сводится к проверке условия непересечения этой гиперплоскости любым из бесконечного множества векторов ai i, i = 1, …, n. Если это условие проверяется эффективно, то алгоритм симплекс-метода также эффективен.

Для того, чтобы понять, насколько эффективен этот измененный для обобщенной задачи линейного программирования симплексный алгоритм, запишем достаточное условие оптимальности базиса в формульном виде. При этом рассмотрим сначала случай, когда нормы pi являются евклидовыми. Именно к такому случаю сводится L-задача оптимального планирования эксперимента.

Пусть B = (a1, …, as) – матрица базиса на текущей итерации, а вектор однозначно находится из системы уравнений Таким образом, относительные разности j для векторов базиса равны нулю (см. шаг 2 алгоритма, приведенного в п. 2.1). Равны они нулю и для любого вектора, конец которого находится на гиперплоскости, проходящей через концы векторов базиса.

Выше этой гиперплоскости (высота измеряется от начала координат) относительные разности отрицательны, а ниже — положительны.

Достаточное условие оптимальности базиса состоит в том, что все векторы из множеств i лежат ниже указанной гиперплоскости. Математически это записывается в виде равенства нулю минимума относительной разности на множестве всех векторов a i, i = 1, …, n:

Из (4.3) получаем 44 4. Задача оптимальной линейной импульсной коррекции и алгоритм ее решения Здесь и максимум реализуется при = Если условие (4.3) не выполняется, т. е. существует такой номер k, что k = max i >1, то в базис вводится любой вектор ak k, для которого не выполняется условие (4.3). Вывод вектора из базиса осуществляется по обычным правилам симплексметода.

Замечание 4.1. Если в рассмотренной задаче коррекции все импульсы одномерные, т. е. si = 1 для всех i, то задача коррекции сводится к задаче линейного программирования, которая возникает при рассмотрении C-задачи планирования эксперимента.

Замечание 4.2. Если полученный указанным алгоритмом оптимальный базис содержит не более одного вектора с ненулевым весом xi из каждого множества i, то этот базис дает решение задачи оптимальной коррекции. Это условие может не выполняться только в случае, когда некоторое множество i содержит плоскость, проходящую через совокупность aij, j = 1, …, ri, векторов из i, имеющих ненулевые веса xij (подробнее об этом см. также в статье [Бахшиян, 1989], где исследуется обобщенная задача линейного программирования). Но в этом случае один из векторов указанной совокупности векторов заменяется некоторой их выпуклой комбинацией — вектором где После этого множество i будет содержать лишь один вектор ai с ненулевым весом, т. е. новый базис дает решение задачи оптимальной коррекции.

Замечание 4.3. Вычисление величины i представляет собой задачу нахождения так называемой двойственной нормы (нормы в сопряженном пространстве) (в нашем случае a = U i, y = i). Для евклидовой нормы двойственная норма также евклидова, а для нормы l1 — это максимальная по модулю компонента вектора a.

Таким образом, для представляющих практический интерес случаев решение оптимальной задачи коррекции осуществляется эффективным алгоритмом, который представляет собой модификацию симплекс-метода, состоящую в особом правиле нахождения вводимого в базис вектора, выбираемого из совокупности множеств i. Поэтому данный метод и получил название «метод генерации столбцов». Этот метод эквивалентен по сложности алгоритму решения обычной задачи линейного программирования, размерность которой (т. е. число линейных уравнений в ограничениях) равна s. При этом приходится оперировать с базисной матрицей размерности ss.

Замечание 4.4. В отличие от обычной задачи линейного программирования, в обобщенной задаче линейного программирования изложенный алгоритм требует, вообще говоря, бесконечного числа итераций (см. рис. 5). Однако за конечное число итераций удается найти приближенное решение с заданной точностью с использованием на каждой итерации оценки для оптимума где max i. Данная оценка является следствием соотношеi ния (2.13).

4.2. ПроеКтная задача идеальной КорреКции Пусть вектор требуемой величины коррекции не задан, а известно лишь, что он лежит в заданной области:

Такой случай имеет место, например, когда при расчетах траектории движения летательного аппарата заранее, еще до запуска, необходимо рассчитать возможный расход топлива. Естественно 46 4. Задача оптимальной линейной импульсной коррекции и алгоритм ее решения поставить задачу нахождения максимального значения суммарного импульса на оптимальную коррекцию при условии, что корректируемый вектор может принимать любые значения из. Таким образом отыскивается Можно показать, что если множество представляет собой выпуклый многогранник, то максимум в (4.5) достигается в одной из вершин этого многогранника.

4.3. КорреКция ПараболичесКой траеКтории летательного аППарата Рассмотрим задачу коррекции траектории летательного аппарата (ЛА), движущегося в гравитационном поле по параболической орбите в некоторой плоскости с известной начальной скоростью. Такое движение может рассматриваться как частный случай движения материальной точки, брошенной под углом к горизонту (см. пример в п. 3.4). Рассмотрим для данного случая задачу идеальной коррекции. Предположим, что коррекция может производиться либо в произвольном направлении, либо в заданном направлении с помощью одного двигателя, жестко связанного с ЛА.

Пусть известны начальная скорость ЛА v0 = v x, v y и вектор корректируемых параметров b.

Требуется решить аналитически задачу коррекции параболической траектории ЛА при изменении двух параметров системы (дальность и высота полета) на величину b.

В момент времени T ЛА должен прийти в требуемую точку, которая является целью полета, однако из-за ошибок исполнения начального импульса реальная траектория ЛА отличается от требуемой и заканчивается в какой-то другой точке. Таким образом, вектор корректируемых параметров b вычисляется как разность параметров для требуемой и реальной траекторий в конечный момент времени.

В случае, когда коррекция проводится в моменты времени ti, i = 1, …, n, задача оптимизации суммарных затрат на коррекцию имеет вид:

4.3. коррекция параболической траектории летательного аппарата Представим ui в виде ui = zii, где zi ui 0 — длина вектора импульса, определяемая нормой, i : i =1 — единичный вектор, определяющий направление импульса.

В результате приходим к следующей задаче:

Тогда задача оптимизации принимает вид Уравнения движения ЛА в произвольный момент времени t записываются в виде Вектор корректируемых параметров определяется как где T — время полета ЛА; x(T ), y(T ) — требуемые значения координат в момент времени T.

Запишем уравнения движения ЛА в момент времени T через значения координат и скорости в произвольный момент времени t :

48 4. Задача оптимальной линейной импульсной коррекции и алгоритм ее решения Рассмотрим два случая, когда импульсы имеют произвольное направление в плоскости движения либо направлены вдоль вектора скорости.

Импульсы направлены произвольно. В данной задаче зависимость r (T) от v (t) нелинейна, поэтому нахождение матрицы влияния сопровождается линеаризацией модели:

где I2 — единичная матрица размерности 22; r (t ) = ( x(t ), y(t )) ;

v(t ) = (v x (t ), v y (t )).

Тогда в дискретном случае а множества i имеют вид т. е. при каждом i концы векторов из множества i лежат на окружности радиуса (T – ti). Если коррекция может проводиться непрерывно в произвольные моменты времени из некоторого интервала [t1, tn], то система множеств i переходит в систему множеств (t ), t1 t tn.

На рис. 6 показан геометрический способ нахождения решения. Несложно заметить, что оптимальным будет приложение одного корректирующего импульса в направлении вектора b в момент времени t = t1.

Найдем величину этого импульса:

4.3. коррекция параболической траектории летательного аппарата Рис. 6. Геометрический способ решения задачи коррекции откуда Импульсы направлены вдоль вектора скорости. В этом случае, как и в предыдущем, корректирующие импульсы ui будут иметь размерность ki = 2. Представим вектор скорости в виде произведения где v(t ) — длина вектора скорости, (t ) = — единичный вектор, определяющий направление вектора скорости.

Из (4.6) следует, что Тогда матрица влияния определяется следующим образом:

50 4. Задача оптимальной линейной импульсной коррекции и алгоритм ее решения Рис. 7. Годограф множеств i (коррекция проводится в заданные дискретные моменты времени) и множеств (t ) (коррекция проводится Рис. 8. Построение решения в задаче оптимальной коррекции:

Рис. 9. Построение оптимального решения в задаче коррекции:

4.3. коррекция параболической траектории летательного аппарата В дискретном случае отсюда следует, что а множества i имеют вид или т. е. каждое из множеств i состоит из двух векторов. На рис. показаны годографы этих множеств в случае, если коррекция может производиться в заданные дискретные моменты времени, и множеств (t ) при t1 t tn, если коррекция может проводиться в любой момент времени.

Возможны следующие случаи расположения вектора b.

А. Вектор b лежит в области годографа множеств i (рис. 8). В этом случае оптимальным решением будет приложение двух корректирующих импульсов (коррекция производится в заданные дискретные моменты времени (см. рис. 8а), оптимальные моменты приложения импульсов — t = t1 и t = t2, импульсы производятся в положительном направлении, т. е. при = +1), либо приложение одного корректирующего импульса, направление которого определяется вектором a(t ) (случай, когда коррекция может проводиться в любой момент времени, см. рис. 8б).

Б. Вектор b лежит вне области годографа множеств i (рис. 9).

В случае, изображенном на рис. 9а, для коррекции ЛА необходимо приложить два корректирующих импульса, один — в направлении вектора a1, т. е. в направлении вектора скорости при t = t1, а второй — в направлении вектора a4, т. е. в момент времени t = t4 в направлении, противоположном направлению вектора скорости.

В случае, изображенном на рис. 9б, для коррекции ЛА также необходимо приложить два корректирующих импульса. Моменты времени и направление импульсов определяются аналогично.

5. гарантирОванные характеристики тОчнОсти 5.1. КритиКа КлассичесКого Подхода К оцениванию точности При довольно естественных условиях, налагаемых на матрицы плана измерений и ковариационную матрицу ошибок, дисперсия оценки Гаусса – Маркова стремится к нулю с увеличением числа измерений n [Эльясберг, 1976]. Однако на практике такой картины не наблюдается. Простой пример приведен в книге [Эльясберг, 1983]. Пусть измеряются показания различных часов одной точности и находится среднее арифметическое всех таких показаний. Такая оценка есть оценка МНК, которая совпадает с оценкой Гаусса – Маркова при отсутствии корреляции между наблюдениями и одинаковой дисперсии. При увеличении числа показаний должно было бы получиться сколь угодно точное время, но на практике этого не произойдет, так как существует корреляция между ошибками в точности хода, обусловленная преемственностью способов изготовления часов, и эта корреляция может меняться в зависимости от времени их изготовления. Поэтому классический подход к вычислению точности оценивания может давать слишком оптимистические результаты.

5.2. гарантирующий Подход К вычислению точности оценивания Этот подход может дать, наоборот, более пессимистические прогнозы точности, но для практики иногда важно «перестраховаться». Поясним сущность этого подхода на рассмотренном выше примере вычисления дисперсии оценки. Пусть известно лишь множество, которому принадлежит ковариационная матрица ошибок измерений. Элементы K этого множества подчиняются условию неотрицательной определенности, которое будем записывать в виде K 0. Тогда в наихудшем с точки зрения величины дисперсии случае получим гарантированную дисперсию Отметим, что гарантированная дисперсия вычисляется при заданном оценивателе x, т. е., например, при заданной весовой матрице, которая уже не может выбираться в соответствии с теоремой Гаусса – Маркова, так как не известна ковариационная матрица ошибок. Вопрос об оптимизации оценивателя рассмотрим ниже. Сейчас выпишем Dmax для простых случаев задания множества. Допустим, что множество определяется следующим образом:

где kij — элементы матрицы K, представляющие собой коэффициенты корреляции между ошибками измерений, число k задано, а дисперсии ошибок kii известны и для простоты записи приняты далее единицами. Тогда где Dk вычисляется без учета условия K 0. Ниже будет показано, что на самом деле это условие при вычислении гарантированной дисперсии можно не учитывать. Тогда или, окончательно, где D0 = xi — дисперсия оценки наименьших квадратов при некоррелированных измерениях (это самый оптимистический 54 5. гарантированные характеристики точности вии, что корреляция может быть произвольной (самый пессимистический случай), т. е. формально при условии kij 1, которое не является ограничением на коэффициенты корреляции, так как выполняется всегда.

Как указывалось в п. 5.1, в неэкзотических случаях D при увеличении числа измерений. Поэтому при большом числе измерений в этих случаях Dmax kD1. Дисперсии D0 и D1 удовлетворяют неравенству Коши – Шварца – Буняковского Из изложенного в разд. 3 материала следует, что, если возможно повторение измерений, то при оптимальном плане измерений здесь будет наблюдаться равенство.

Отметим более общий случай [Бахшиян и др., 1980], когда коэффициенты корреляции принадлежат многомерному прямоугольному параллелепипеду, т. е. элементы kij удовлетворяют неравенствам где kij — элементы некоторой известной номинальной корреляционной матрицы K, а неотрицательные элементы wij задают разброс около K (kii = 1, wii = 0 при всех i). В этом случае где вектор x+ состоит из компонентов xi. Критерий достижимости равенства Dmax = Dw приведен в книге [Бахшиян и др., 1980], но этот критерий является трудно проверяемым. Однако имеется простое достаточное условие достижимости равенства в последнем соотношении, которое имеет вид при некоторой диагональной матрице. В случае K = = I, где I — единичная матрица, получаем W + I 0 и Dmax = Dw.

Отсюда следует упомянутое выше равенство Dmax = Dk.

5.3. сравнение решений задач оптимального оценивания в двух простейших случаях… 5.3. сравнение решений задач оПтиМального оценивания в двух Простейших случаях При гарантирующеМ и КлассичесКоМ Подходах Покажем удивительное совпадение состава оптимальных измерений при двух различных подходах к модели ошибок. Рассмотрим сначала классический подход.

5.3.1. оптиМизация гарантированной диСперСии D Задача оптимизации алгоритма оценивания с точки зрения минимизации гарантированной дисперсии имеет вид Выбирая различные скалярные параметры lj, j =1,, k m, находим для них различные оптимальные векторы-оцениватели. Составляя из этих векторов матрицу оценивателя X размерности km, можно найти, согласно изложенному ранее, множество весовых матриц таких, что каждой матрице соответствует МНК, дающий ту же оценку, что и оцениватель X.

Рассмотрим задачу минимизации D1 (напомним, что D0 соответствует единичной корреляционной матрице и найдена в п. 5.2). Имеем, согласно (5.2), и величина D1 связана с величиной D0 соотношением (см. п. 5.2).

Отсюда следует важный методологический вывод [Бахшиян, 1970]. В случае повторяющихся некоррелированных измерений имеется m оптимальных точек повторения измерения, которые совпадают с оптимальными точками проведения измерений ( xi 0) для задачи минимизации гарантирующей дисперсии D1, полученной при возможности произвольной корреляции между измерениями. При этом D1 = nD0, а число повторений измерений ri пропорционально модулю коэффициента xi, соответствующего i-му измерению в случае произвольной корреляции (см. п. 3.3).

56 5. гарантированные характеристики точности 5.3.2. МиниМакСная задача оценивания при ограниченных по Модулю ошибках изМерений Такая задача рассматривалась М. Л. Лидовым [Лидов, 1964] и названа им «схемой бортиков». Предполагается, что ошибки измерений ограничены по модулю:

где Mi — известные числа. Тогда гарантированная ошибка измерений равна а ее минимизация сводится к задаче линейного программирования того же типа, что и задачи в п. 3.3, 5.3.1:

Рассмотрим частный случай, когда Mi = M при всех i. Этот случай можно трактовать как допущение о пропорциональности максимально возможных ошибок измерений их дисперсиям, которые в п. 5.3.1 приравнены к единицам. Тогда задачи линейного программирования, решаемые при нахождении величин l и 1 (см. п. 5.3.1), дают один и тот же оптимальный оцениватель, при этом Отметим, что вместо условий на ошибки измерений можно было бы рассматривать аналогичные условия на их математические ожидания:

При этом гарантированное математическое ожидание равно а задача минимизации Emax по линейному несмещенному оценивателю приводит к той же по виду задаче линейного программирования, что и задача нахождения l (см. выше).

5.4. вычисление гарантированных характеристик точности оценивания… Замечание 5.1. Можно показать, что решение указанной задачи линейного программирования дает оптимальный алгоритм оценивания в классе всех (как линейных, так и нелинейных) несмещенных алгоритмов в «схеме бортиков» [Матасов, 1988а, б]. Это позволяет более прагматично отнестись к казалось бы экзотической «схеме бортиков» и философски посмотреть на проблему использования всей информации в том случае, когда информация поставляется «игроком», сопротивляющимся получению хорошей оценки.

5.4. вычисление гарантированных хараКтеристиК точности оценивания При наличии неМоделируеМых возМущений 5.4.1. опиСание Моделей движения и изМерений Пусть X = X (t ) — вектор, полностью определяющий траекторию движения системы в заданной системе координат.

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

Пусть изменение вектора X на интервале [0, T] описывается системой дифференциальных уравнений где u = u(t ) — вектор возмущений, значение которого мало, но не известно.

Возмущение u(t) далее будем называть немоделируемым.

Ввиду неопределенности этого возмущения при решении задачи определения движения системы могут быть использованы лишь модельные уравнения движения, которые получаются из (5.3) при нулевом немоделируемом возмущении:

где 0r — нулевой вектор. Пусть 58 5. гарантированные характеристики точности — начальные условия для некоторого заданного момента времени t0, а — решение системы (5.4) при этих начальных условиях. Задача оценивания траектории по результатам измерений сводится к нахождению оценки X 0 вектора X0, после чего оценка вектора состояния определяется согласно (5.6):

Если скалярный параметр траектории l есть некоторый функционал от X(t), т. е.

то оценка скалярного параметра (5.8) находится при известной оценке (5.7) по формуле В частности, параметр l может быть функцией состояния в некоторый заданный момент времени :

Подставляя (5.6) в (5.8) или (5.10), получаем в результате зависимость вида Поэтому соотношение (5.9) для оценки параметра l эквивалентно соотношению, получаемому при использовании функции (5.11):

Перейдем теперь к описанию модели измерений. Пусть на некотором множестве [0, T ] производятся измерения Здесь y(t) — скалярная функция измеренных значений; (·) — функция известного вида; y(t) — ошибка измерений. Множество на практике состоит из конечного числа моментов времени ti, которые будем полагать упорядоченными:

5.4. вычисление гарантированных характеристик точности оценивания… Иногда нам будет удобно рассматривать общий случай, когда это множество может включать произвольные числовые множества из [0, T]. Если эти множества — интервалы, то будем говорить о непрерывных измерениях, а в случае (5.14) — о дискретных измерениях.

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

5.4.2. линеаризация Модельных уравнений Далее будем полагать, что уравнения (5.3), (5.13) могут быть линеаризированы около некоторой опорной траектории, характеризуемой фазовым вектором X (t ) Наилучшей была бы линеаризация около истинной траектории движения, однако она не известна. Поэтому опорная траектория уточняется путем ее определения по результатам измерений. В итоге ошибка линеаризации будет иметь порядок величины X (t ) - X (t ), где — евклидова норма. Этой ошибкой линеаризации можно либо пренебречь, либо при линеаризации уравнений (5.3), (5.13) учесть эту ошибку путем соответствующего добавления ее к немоделируемым возмущениям и ошибкам измерений. При линеаризации выражения (5.10) для оцениваемого параметра добавляется ошибка линеаризации l.

В результате вместо соотношений (5.3), (5.13), (5.10) получаем для вектора линейные уравнения движения и измерений:

и выражение для оцениваемого параметра:

60 5. гарантированные характеристики точности (l = 0, если ошибка линеаризации пренебрежимо мала). Здесь R(t), F(t) — известные матричные функции размерности mm и mr соответственно; h(t ), g (t ) — известные векторm функции; a — известный вектор. Функции R(t) и F(t) вычисляются по формулам где функция f(·) определена в (5.3).

Представим модель измерений (5.18) как сумму двух слагаемых: первое зависит лишь от вектора состояния в момент t0, второе есть суммарная ошибка в модели (5.18), возникающая вследствие ошибок измерения и немоделируемых возмущений.

Для этого выпишем решение системы (93) при начальном условии q(t0) =. Имеем согласно формуле Коши:

где Q(t) — фундаментальная матрица размерности mm, определяемая из системы дифференциальных уравнений Im — единичная матрица размерности mm, m — размерность вектора состояния.

5.4. вычисление гарантированных характеристик точности оценивания… Подставляя (5.22) в (5.18), получаем новую запись для модели измерений и оцениваемого параметра:

где (t) — суммарная ошибка исходных данных, обусловленная ошибками измерений и немоделируемыми возмущениями, — аналогичная ошибка в модели оцениваемого параметра.

Для рассматриваемого случая дискретных измерений запишем модель (5.23) в векторном виде:

где Y — вектор измерений с компонентами y(ti), i = 1, …, n, — матрица размерности mn, которую будем называть матрицей плана измерений.

5.4.3. вычиСление гарантированной ошибки линейного оценивания Будем рассматривать линейную оценку вида где оцениватель xi, i = 1, …, n, удовлетворяет условию несмещенности Ошибка оценивания параметра l для несмещенного оценивателя имеет вид 62 5. гарантированные характеристики точности где Далее будем пренебрегать ошибкой линеаризации l и считать, что g (t ) 0. При необходимости эти слагаемые могут быть учтены. Для того, чтобы правильно учесть влияние u(t), нужно собрать множители при u(t), соответствующие одному моменту времени.

Это наиболее просто делается в случае > tn, для которого Здесь В общем случае формулы записываются в виде где (x) — функция Хевисайда:

Замечание 5.2. Вектор Yi(t) зависит от t только на отрезке [tk – 1, tk], которому принадлежит момент.

Замечание 5.3. Можно показать, что в рамках используемого линейного приближения равенство не зависит от выбора начального момента t0.

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

Далее будем полагать l = 0. Пусть lu max — максимальная ошибка прогноза, обусловленная погрешностью u(t) при условии, что эта погрешность принадлежит некоторому известному множеству Mu:

5.4.4. вычиСление lu max для некоторых Mu 1. Пусть известно лишь, что u(t) — неслучайная величина, ограниченная по евклидовой норме:

где (t) — некоторая известная функция. Такое ограничение естественно принять при допущении о том, что немоделируемое возмущение в каждый момент времени обусловлено, например, истечением газов в космическом аппарате в неизвестном направлении.

Используя те же идеи, что и в разд. 4, получаем Гарантированная ошибка получается при 2. Пусть известно, что компоненты вектора u(t) ограничены по модулю:

64 5. гарантированные характеристики точности Такое ограничение естественно принять при допущении о том, что немоделируемое возмущение в каждый момент времени обусловлено, например, истечением газов в космическом аппарате в заданных направлениях. Тогда получаем где Z i (t )+ — вектор, компоненты которого состоят из модулей компонент вектора Zi(t); (t) — вектор с компонентами (t)j.

5.5. вычисление гарантированных хараКтеристиК точности суММарной ошибКи оценивания В соответствии с (5.31) ошибка оценивания параметра l при отсутствии ошибок линеаризации равна где ly — ошибка, обусловленная погрешностями измерений, которая вычисляется по формуле (5.32), а ошибка lu возникает вследствие неучета немоделируемых возмущений и описывается формулой (5.38). Допустим, что математические ожидания ошибок измерений y(ti) и их коэффициенты корреляции ограничены по модулю, т. е.

значение дисперсии ошибки l равно Здесь lu max вычисляется в зависимости от принятых допущений по одной из формул п. 5.4.4, величина Dk находится из (5.2).

Если предположить, что ошибки измерений ограничены по модулю, т. е. y (ti ) M i, то из (5.44) получаем максимальное значение ошибки оценивания:

5.6. оптимизация гарантированных характеристик точности методом генерации столбцов Здесь lu max как нелинейная функция от коэффициентов xi оценивателя определена выражением (5.41) или (5.43) в зависимости от принятых допущений. Далее будет показано, как может быть решена задача нахождения оптимального несмещенного оценивателя путем минимизации выражения (5.48).

5.6. оПтиМизация гарантированных хараКтеристиК точности МетодоМ генерации столбцов Рассмотренные выше оптимальные задачи оценивания сводятся либо к задачам линейного программирования, число ограничений которых равно размерности m вектора оцениваемых параметров, либо к задачам обобщенного линейного программирования с числом ограничений, пропорциональным m. Трудоемкость таких задач пропорциональна числу их ограничений.

Более сложные практические задачи (например, задача минимизации выражения (5.48) по xi при наличии условия несмещенности) не могут быть сведены к таким эффективно решаемым при небольших m задачам.

Рассмотрим общий подход, позволяющий решать достаточно широкий класс задач методом генерации столбцов, упомянутым в разд. 4. Недостатком этого подхода является то, что число ограничений решаемой обобщенной задачи линейного программирования равно n + 1, где n — число измерений. Это приводит к большой трудоемкости метода генерации столбцов при больших n.

Рассмотрим сначала задачу минимизации гарантированной ошибки оценивания в более общей постановке, чем ранее.

Пусть задано симметричное относительно нуля выпуклое замкнутое ограниченное множество M возможных значений вектора ошибок измерений. Рассмотрим минимаксную задачу оценивания которая дает минимальное значение ошибки в наименее благоприятном случае. Такой вид имеет, например, задача о минигарантированные характеристики точности мизации гарантированной ошибки (5.48). Задачу (5.49) можно переписать также в виде где Будем при этом предполагать, что функция f(x), которая равна гарантированной ошибке оценивания при заданном оценивателе x, может быть вычислена либо аналитически, либо достаточно эффективно численными методами.

Переставим операции минимума и максимума (это возможно согласно теореме фон Неймана) и учтем, что согласно теории двойственности Тогда получим Покажем, что задачу нахождения l удается решить эффективно методом генерации столбцов для обобщенного линейного программирования, изложенным в разд. 4. Для этого в соответствии с теоремой Каратеодори представим каждый элемент выпуклого множества M как выпуклую комбинацию неизj) вестных n + 1 точек M, j = 1, …, n + 1. Тогда выписанная выше задача максимизации по запишется в виде при условиях где Hi — n-мерные столбцы матрицы H.

Это есть обобщенная задача линейного программирования, размерность которой (т. е. число уравнений в ограничениях) равна n + 1. Неограниченные по знаку переменные i можно всегда принимать базисными при решении обобщенной задачи линейного программирования методом генерации столбцов. Это следует из того, что, если соответствующие столбцы попадают 5.6. оптимизация гарантированных характеристик точности методом генерации столбцов в базис, то они оттуда не убираются при использовании метода генерации столбцов (см. (2.9)).

( — скаляр), определяемый из условия (см. разд. 2):

где B — базисная матрица размерности (n + 1)(n + 1), соответствующая текущему базису. Тогда достаточным условием оптимальности текущего базиса является выполнение равенства вида (4.3):

и в базис вводится вектор условий (-,1), для которого >.

Таким образом, для проверки условия оптимальности и введения в базис вектора условий нужно уметь решать подзадачу которая по сделанному выше предположению о функции f(x) решается эффективно (например, аналитически). Поэтому вычисление величины l методом генерации столбцов для задачи (5.52) эффективно. К тому же, согласно (2.14), имеется оценка близости текущего значения b целевой функции к оптимуму:

Можно показать, что скаляр равен оптимальному значению l, а в качестве оптимального оценивателя x можно взять вектор, удовлетворяющий достаточным условиям оптимальности. Это следует из того, что двойственный вектор к двойственной задаче (5.52) является вектором оптимизации в прямой задаче (5.49).

литература [Бахшиян, 1970] Бахшиян Б. Ц. Выбор оптимальных моментов независимых траекторных измерений // Космические исследования. 1970. Т. 8.

[Бахшиян, 1983] Бахшиян Б. Ц. Представление весовых матриц, определяющих заданную оценку наименьших квадратов // Навигационная привязка и статистическая обработка космической информации. М.:

Наука, 1983.

[Бахшиян, 1989] Бахшиян Б. Ц. Критерии оптимальности и алгоритмы решения вырожденной и обобщенной задач линейного программирования // Экономика и математические методы. 1989. Т. 28. № 2. С. 314–324.

[Бахшиян, Соловьев, 1998] Бахшиян Б. Ц., Соловьев В. Н. Теория и алгоритмы решения задач L- и MV-оптимального планирования эксперимента // Автоматика и телемеханика. 1998. № 8. С. 80–96.

Бахшиян и др., 1980] Бахшиян Б. Ц., Назиров Р. Р., Эльясберг П. Е. Определение и коррекция движения. М.: Наука, 1980.

[Бахшиян и др., 2000] Бахшиян Б. Ц., Матасов А. И., Федяев К. С. О решении вырожденных задач линейного программирования // Автоматика и телемеханика. 2000. № 1. С. 105–117.

[Войсковский, 1998] Войсковский М. И. Симплексный алгоритм решения задачи MV-оптимального планирования эксперимента: Препринт ИКИ РАН. М.: ИКИ РАН, 1998. Пр-1979.

[Данциг, 1966] Данциг Дж. Линейное программирование, его применения и обобщения. М.: Прогресс, 1966.

[Лидов, 1964] Лидов М. Л. К априорным оценкам точности определения параметров по методу наименьших квадратов // Космич. исслед. 1964.

Т. 2. № 5. С. 713–718.

[Лидов, 1971] Лидов М. Л. Математическая аналогия между некоторыми оптимальными задачами коррекции траекторий и выбора состава измерений и алгоритмы их решения // Космич. исслед. 1971. Т. 9. № 5.

С. 687–706.

[Лэсдон, 1975] Лэсдон Л. С. Оптимизация больших систем. М.: Наука, 1975.

[Матасов, 1988а] Матасов А. И. Об оптимальности линейных алгоритмов гарантирующего оценивания. Часть I // Космич. ислед. 1988. Т. 26.

№ 5. С. 643–653.

[Матасов, 1988б] Матасов А. И. Об оптимальности линейных алгоритмов гарантирующего оценивания. Часть II // Космич. ислед. 1988. Т. 26.

№ 6. С. 807–812.

[Эльясберг, 1976] Эльясберг П. Е. Определение движения. М.: Наука, 1976.

[Эльясберг, 1983] Эльясберг П. Е. Измерительная информация: сколько ее нужно? как ее обрабатывать? М.: Наука, 1983.

[Elfving, 1952] Elfving G. Optimum Allocation in Linear Programming // Annals of Mathematical Statistics. 1952. V. 23. P. 255.

[Rao, 1975] Rao C. R. On a Unified Theory of Estimation in Linear Models — a Review of Recent Results // Perspectives in Probability and Statistics. L.

etc.: Acad. Press, 1975. P. 89–104.

доПолнительная литература Бажинов И. К., Почукаев В. Н. Оптимальное планирование навигационных измерений в космическом полете. М.: Машиностроение, 1976.

Белоусов Л. Ю. Определение оптимальных моментов измерений // Космич.

исслед. 1969. Т. 7. № 1. С. 28–34.

Белоусов Л. Ю., Комаров В. А. Некоторые общие результаты в задаче определения оптимальных моментов измерения для различных моделей ошибок измерений // Космич. исслед. 1970. Т. 8. № 3. С. 452–453.

Гольштейн Е. Г. Выпуклое программирование. Элементы теории. М.: Наука, 1989.

Войсковский М. И. Симплексный алгоритм поиска E-оптимальных планов // Изв. РАН. Теория и системы управления. 2001. № 2. C. 70–74.

Войсковский М. И., Меринов И. Е. Симплексный алгоритм решения минимаксной задачи оценивания: Препринт ИКИ АН СССР. М.: ИКИ АН СССР, 1990. Пр-1697.

Демиденко Е. З. Линейная и нелинейная регрессии. М:. Финансы и статистика, 1981.

Ермаков С. М., Жиглявский А. А. Математическая теория оптимального эксперимента. М.: Наука, 1987.

Иоффе А. Д., Тихомиров В. М. Теория экстремальных задач. М.: Наука, 1974.

Лидов М. Л. Минимаксные методы оценивания. М.: Изд-во ИПМ РАН, Лидов М. И., Бахшиян Б. Ц., Матасов А. И. Об одном направлении в проблеме гарантирующего оценивания (обзор) // Космич. исслед. 1991.

Т. 29. № 5. С. 659–684.

Лидов М. И., Матасов А. И. Об одном обобщении задачи о «наихудшей корреляции» // Космич. исслед. 1989. Т. 27. № 3. С. 454–456.

Матасов А. И. Задачи гарантирующего оценивания: Курс лекций. М.: Издво МГУ, 2009.

Мину М. Математическое программирование. Теория и алгоритмы. М.:

Наука, 1989.

Муртаф Б. Современное линейное программирование. М.: Мир, 1984.

Федоров В. В. Теория оптимального эксперимента. М.: Наука, 1971.

Pukelsheim F. Optimal design of experiments. N. Y.: J. Wiley and Sons, 1993.

1. Задача оценивания параметров движущейся системы 1.2. Несмещенный алгоритм оценивания.................... 1.5. Линейный несмещенный алгоритм 1.6. Вычисление точности оценок при известной ковариационной матрице ошибок измерений.

1.7. Эквивалентность множеств всевозможных линейных несмещенных оценок при различном выборе мешающих параметров........... Основные сведения из теории линейного программирования..... 2.1. Постановка задачи и симплекс-метод................... 2.2. Оценка близости текущего решения к оптимальному 2.3. Построение допустимого базисного решения............ 2.4. Геометрический способ решения двумерной задачи линейного программирования специального вида........ 3. Классические задачи планирования эксперимента и их сведение к задачам оптимальной идеальной импульсной коррекции 3.3. Скалярная задача планирования эксперимента и алгоритм 3.4. Оценивание параметров параболической траектории 3.5. Сведение L-задачи к задаче оптимальной линейной 3.6. Сведение MVs-задачи к параметрической задаче оптимальной линейной импульсной коррекции 3.6.1. Сведение MVs-задачи к проектной задаче 3.6.2. Редукция задачи (3.22) к задаче многомерной 4. Задача оптимальной линейной импульсной коррекции 4.1. Задача оптимальной коррекции 4.2. Проектная задача идеальной коррекции................. 4.3. Коррекция параболической траектории летательного 5. Гарантированные характеристики точности................... 5.1. Критика классического подхода к оцениванию 5.2. Гарантирующий подход к вычислению точности 5.3. Сравнение решений задач оптимального оценивания в двух простейших случаях при гарантирующем 5.3.1. Оптимизация гарантированной дисперсии D1....... 5.3.2. Минимаксная задача оценивания при 5.4. Вычисление гарантированных характеристик точности оценивания при наличии немоделируемых возмущений... 5.4.1. Описание моделей движения и измерений.......... 5.4.2. Линеаризация модельных уравнений.............. 5.4.3. Вычисление гарантированной ошибки линейного 5.5. Вычисление гарантированных характеристик точности 5.6. Оптимизация гарантированных характеристик точности