Оглавление

# Поиск оптимальных решений обратных задач

## Основы UNEX

Основной метод решений обратных задач в UNEX - это метод наименьших квадратов (МНК), который позволяет минимизировать функционал, оценить ошибки полученных параметров, степень корреляции параметров и др. Такие задачи в UNEX выполняются с помощью команды MINIMIZE и ее разновидностей (например ROBUSTM). Но обычно мы имеем дело с нелинейным МНК и поэтому получаемое нами решение зависит от начального приближения, т.е. мы сталкиваемся с понятиями локальных и глобального минимумов функционала. Начав минимизацию, будучи в одном минимуме, мы не сможем прийти в другой, возможно более глубокий минимум (для примера см. рис. 2). Поэтому для поиска начального приближения (и не только) полезно пользоваться другим семейством команд - SEARCH. Сначала будут описаны основные особенности команд, не затрагивая конкретных методов исследования молекул, а потом будут даны описания команд, относящихся к конкретным методам и способам исследования молекул.

## MINIMIZE

Эта команда минимизирует заданный функционал. Метод минимизации - МНК или метод золотого сечения или же комбинированный метод (и то и другое). Включить тот или иной метод можно ключем MinMethod в поле BASE. По умолчанию используется комбинированный метод, заключающийся в том что основную работу выполняет МНК, а в случае плохой обусловленности задачи - подключается метод золотого сечения как вспомогательный. Такая связка методов совершенно безотказна (конечно если модель задана правильно), хотя и довольно медленная при работе оптимизации методом золотого сечения.
В качестве математической основы МНК применяется метод Левенберга-Макгварда [4]. Линеаризованная проблема решается методом SVD.
В начале МНК-минимизации проводится локальный анализ чувствительности первого рода.Выводится таблица (Local first-order sensibility analysis) с двумя столбцами: номера групп и относительные коэффициенты чувствительности. Эти коэффициенты показывают относительную степень влияния каждой группы параметров на модель в НАЧАЛЬНЫЙ (локальный анализ!) момент минимизации. По этим числам можно судить о том, насколько далеко находится начальное значение параметра (или группы параметров) от оптимального и насколько сильно этот параметр (или группа) будет влиять на изменение функционала.
Во время минимизации на экран и в лог выводится в общем случае три столбца чисел:

1. Текущей номер МНК-итерации.
2. Текущее значение электронографического функционала относительно его начального значения.
3. Добавочный коэффициент в методе Левенберга.

Если в минимизации учавствуют электронографические данные, то выводится также текущий общий Rf-фактор.
По динамике изменения коэффициента Левенберга можно проследить устойчивость нахождения минимума: если этот коэффициент постоянно снижается, то система устойчиво движется к минимуму. Если значение коэффициента меняется в разные стороны, то наблюдаются осцилляции и система скорее всего идет по пологому оврагу или уже находится в окресности пологого минимума. Сами по себе такие осцилляции не страшны, метод Лвенберга все равно подойдет максимально близко к минимуму, но в тоже время нет особого смысла тратить много итераций в таких ситуациях - лучше пробовать изменить составы групп и включить/выключить варьирование других параметров чтобы модель была более устойчивой.
Если слишком быстро начал увеличиваться добавочный коэффициент Левенберга, то задача плохо обусловлена и скорее всего система находится в пологом месте на гиперповерхности. Если метод Левенберга-Макгварда в пологом минимуме исчерпает свои возможности, то (если MinMethod не установлен в другое щначение) производится попытка подойти ближе к минимуму с помощью покоординатного спуска методом золотого сечения. После отработки метода золотого сечения, в случае когда число допустимых итераций (ключ MaxIter в поле BASE) еще не исчерпано, продолжается минимизация методом Левенберга-Макгварда. Обычно это свидетельствует о математической "слабости, вырожденности, плохой обусловленности" модели, и лучше перераспределить параметры в группах. Ускорить (или наоборот, замедлить для устойчивости) ход минимизации можно с помощью параметра damp в поле BASE. Все МНК-добавки к оцениваемым параметрам умножаются на это число. По умолчанию это число изменяется сигмоидально в ходе минимизации. На начальных итерациях damp-фактор горазо меньше еденицы (для придания устойчивости ходу минимизации). На последней итерации damp равен еденице. Если вы уверены что задача устойчива, то можно зафиксировать damp на одном значении, например damp=1.0, чтобы быстро прошла минимизация. Есть вариант damp=linear, когда damp-дамп фактор линейно увеличивается в ходе минимизации.
Минимизация останавливается в трех случаях:

1. Допустимое количесвто итераций выполнено. В этом случае рекомендуется увеличить допустимое количество итераций ключем MaxIter в поле BASE, или позаботиться об "усилении" модели, перераспределив группы параметров.
2. Найден минимум функционала, найдены оптимальные параметры. Это найлучшый вариант.
3. Относительное значение функционала пренебрежимо мало, найдены оптимальные параметры.

Второй вариант достижим при одновременном выполнении условий:

1. Максимальная по модулю производная параметра по функционалу меньше критического значения.
2. Максимальное по модулю изменение производной меньше критического значения.
3. Изменение целевого функционала меньше критического значения.
4. Максимальная итерационная добавка меньше критич. значения.
5. Максимальное изменение добавки меньше критич. значения.

Критические значения, принятые по умолчанию (их можно менять с помошью специальных ключей в поле BASE) для вышеуказанных параметров в том же порядке:

    1. 0.000100
    2. 0.000010
    3. 0.000010
    4. 0.000001
    5. 0.000001

В конце минимизации выводится некоторая статистическая информация, матрица коэффициентов корреляции, таблица оцениваемых параметров и эллипсоид (в общем случае - гиперэллипсоид) ошибок. В состав стат. информации входит:

1. абсолютное значение функционала - X^2,
2. конечный параметр Левенберга-Макгварда - La,
3. параметр Goodnes-of-Fit (имеет значение, если заданы адекватные ошибки измерения экспериментальных данных) - GOF,
4. Если минимизация шла по ГЭ данным: структурные R-факторы по каждой кривой sM(s) в отдельности, полный структурный R-фактор.
5. Степень обусловленности линеаризованной системы уравнений - отношение найбольшего собственного числа к найменьшему.

В таблице параметров представлены следующие величины: номер группы/тип параметра/старое значение первого параметра в группе/новое значение первого параметра в группе/утроенное значение среднеквадратичного отклонения/относительная ошибка в процентах/значение частной производной функционала по данной переменной/частный коэффициент корреляции для данной группы. Вывод значений углов и их ошибок подчиняется ключу AngleUnits. По умолчанию углы (и их ошибки) выводятся в градусах.
Если в матрице корреляций есть хотя бы одно число, превышающее по модулю 0.5, то будет выведена дополнительная таблица коэффициентов коррелиции. В этой таблице коэффициенты корреляций разбиты на группы: более 0.5, более 0.6, более 0.7, более 0.8 и более 0.9. Такая таблица удобна для быстрого обзора больших коэффициентов корреляции.
Эллипсоид ошибок выводится в виде набора ортонормированных векторов (соответствующим полуосям эллипсоида) и длин полуосей. Смысл эллипсоида: решение МНК-задачи находится в пределах эллипсоида с вероятностью 99.73%.
Вот пример вывода гиперэллипсоида (пятимерного):

     X^2 hyperellipsoid (axes vectors) for P=99.73%
     ------------------------------------------------------------------------------
        NN   |    Length    |  Vector components --->
     ------------------------------------------------------------------------------
         1    2.18685e-003  : -0.002204  0.004154  0.998881 -0.007848  0.046394
         2    4.67001e-003  :  0.061441  0.001608  0.046089 -0.039742 -0.996252
         3    1.16254e-002  : -0.061122  0.000128  0.009701  0.997152 -0.043098
         4    1.83963e+000  :  0.991919  0.092900 -0.000430  0.063334  0.058777
         5    4.80662e+000  :  0.092632 -0.995666  0.004203  0.005940  0.004063
     ------------------------------------------------------------------------------

В этой таблице указаны номер каждого вектора (вектор соответствует оси эллипсоида), длина вектора (что соответствует длине полуоси эллипсоида) и компоненты вектора в системе координат варьировавшихся параметров.
В сложных случаях, когда количество оцениваемых параметров велико и/или наблюдаются значительные коррелиции между параметрами, имеет смысл варьировать не сами параметры, а их ортогональные линейные комбинации со специально подобранными коэффициентами. Это достигается путем включения в поле BASE ключа MinOrthoParams=1. Варьирование ортогональных параметров численно более устойчиво. В обычных ситуациях структурного анализа, при использовании дампирования и стандартных критериев остановки минимизации функционала, варьирование ортогональных параметров приводит к тем же результатам что и варьирование обычных параметров.

## ROBUSTM

Это команда для проведения минимизации, устойчивой по отношению к грубым промахам в экспериментальных данных. На практике это достигается с помощью итеративного изменения весов экспериментальных данных по схеме "Tukey's bisquare weights". Команда ROBUSTM имеет точно такой же синтаксис как и у команды MINIMIZE.

## SEARCH

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

    SEARCH=способ_поиска,тип_функционала,....

Способы поиска:

1. SCAN - простой перебор значений параметров в заданных пределах с заданным шагом (далее - сканирование).
2. RAND - метод Монте-Карло (еще говорят "рандомизация"). Параметры выбираются случайным образом в заданных пределах.

Типы функционала все те же, что и в команде MINIMIZE.
Далее будут даны описания различных способов работы команды SEARCH:
SCAN
Это простой способ тестирования значений функционала при различных параметрах. Параметры, их диапазоны изменения и количества точек указываются в специальном поле, где находится информация для сканирования:

    SEARCH=SCAN,тип_функционала,<scan_tag>,</scan_tag>,....

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

    <scan_tag>
    1   1.00    5     2.00
    2   1.57   10     1.87
    5  80.00   50   115.00
    3   0.01   10     0.10
    4   0.02   10     0.10
    </scan_tag>

В каждой строчке указано: номер группы, найменьшее значение параметра, количество точек на этот параметр, найбольшее значение параметра. Здесь приведен пример задания на пятимерное сканирование. Значения сканируемых параметров будут присваиваться первым представителям указанных групп, остальные члены группы будут иметь значения отличающиеся от значения первого члена на начальную разницу, т.е. разности параметров внутри одной группы при сканировании сохраняются. Число размерностей не ограничено. Число точек на каждую размерность не ограничено.
Во время сканирования на консоль выводится текущий процент выполненной работы и скорость сканирования. Убрать эту информацию можно ключем ShowScanInfo в поле BASE. По окончании сканирования в лог выводится (если PrintScanResults в поле BASE не установлено в "no") информация: таблица полученных значений функционала (или R-факторов), и таблица с параметрами отвечающими минимальному значению функционала. При этом таблица выглядит следующим образом. Если задается одномерное сканирование, то в лог выводится два столбца чисел - значние параметра и значение функционала (или R-фактора). Если используется двумерное сканирование - то в лог выдается матрица значний функционала (или R-фактора). Если же используется сканирование большей размерности - то выдаются матрицы (называемые блоками), которые соответствуют сечениям многомерной поверхности. Пример результатов двумерного сканирования в графическом виде представлен на рис. 1.
Рис. 1: Поверхность значений функционала в зависимости от двух параметров.
RAND
Метод Монте-Карло. Значения переменных выбираются случайным образом и вычисляется значение функционала. Такая процедура происходит многократно. Останавливается эта процедура по истечении времени указанном с помощью ключа SrchTime в поле BASE. Выдаются в лог значения параметров, отвечающих найменьшему значению функционала. При RAND поиске, в отличии от варианта SCAN, таблица значений функционала при разных значениях параметров не выводится.
Особенности команды SEARCH:

1. Время работы команды SEARCH контролируется ключем SrchTime в поле BASE. SrchTime - единственный способ контроля длительности работы команды SEARCH=RAND.
2. Если среди варьируемых командой SEARCH параметров есть больше одной концентрации (в долях еденицы) и их диапазоны, например от 0.0 до 1.0, то ясно что невозможны комбинации, когда обе концентрации больше 0.5. В этом случае UNEX отнормирует концентрации так, чтобы сумма концентраций всех молекул в модели не превышала еденицы. По этой же причине, таблица значений функционала (или R-фактора), выводимая после команды SEARCH=SCAN может иметь искаженный вид. Однако определенные таким образом оптимальные параметры все равно будут корректными.
3. После того как отработала команда SEARCH, модель имеет найлучшие (найденные SEARCH) параметры.

# Обратная задача газовой электронографии

Обычно в структурном анализе газовой электронографии минимизация функционала осуществляется методом найменьших квадратов (МНК). Но часто сложность задачи такова, что одним методом МНК обойтись нельзя. Например, если заранее не известно какой конформер лучше всего описывает эксперимент, то приходится перебирать "вручную" возможные решения. Особенно ощутим такой перебор при конформационных пространствах большой размерности. Здесь в значительной мере помогает процедура сканирования структурного R-фактора по различным параметрам.

## MINIMIZE

На уровне использования экспериментальной электронографической информации, в программе реализована "классическая" схема, когда используются экспериментальные кривые sM(s). Так что слово GEDINT в команде

    MINIMIZE=GEDINT,...

означает минимизацию на кривых sM(s), т.е. минимизацию функционала:
Где: i - номер кривой sM(s), j - номер точки на i-ой кривой sM(s), Scale - масштабный множитель, w - вес (равен обратному квадрату ошибки), символы "e" и "t" обозначают эксперимент и теорию соответственно.
Перед тем как начать минимизацию, вам нужно определится с тем, какие параметры вы хотите варьировать:

1. Масштабные множители Scale. По умолчанию (и если в команде MINIMIZE явно не указывать номеров групп) они всегда варьируются. При желании можно выключить их варьирование просто поставив у соответствующей кривой VarSc=0 (см. секцию про ввод интенсивностей).
2. Мольные доли компонентов в долях еденицы. Номера групп таких параметров задаются в собственных полях соответствующих молекул. Помните, что в смеси из N молекул, можно варьировать не больше N-1 концентраций.
3. Геометрические параметры. Номера групп геометрических параметров задаются во второй части (часть описания переменных) z-матриц молекул.
4. Среднеквадратичные амплитуды колебаний. Номера групп задаются в полях ввода амплитуд.
5. Параметры потенциальной функции. По умолчанию не варьируются. При вводе потенциала можно задавать номера групп параметрам потенциальной функции.

Основные правила объединения параметров в группы:

1. Нельзя группировать разнородные параметры. Т.е. нельзя в составе одной группы варьировать:
   - связи и углы,
   - связи и амплитуды,
   - углы и амплитуды,
   - концентрации и связи/углы/амплитуды,
   - коэффициенты потенциальной функции и связи/углы/амплитуды.
2. Во избежание сильных корреляций и как следствие больших ошибок и смещенных оценок, не рекомендуется варьировать близкие длины связей (меньше 0.03-0.05 ангстрема) в разных группах. Вместе с этим, не рекомндуется объединять в одну группу связи, вносящие существенно разный вклад в дифракционную картину, например C-H и C=O. В этом случае C-H лучше всего зафиксировать. Аналогично с валентными углами: близкие по величине углы, построенные на одинаковых (или близких) связях, контролируют близкие несвязанные межъядерные расстояния. Такие углы крайне не рекомендуется варьировать в разных группах. Например, нехорошо варьировать валентные углы в бензольном кольце в разных группах.
3. Амплитуды пар атомов, расстояния между которыми близки, нужно варьировать в одной группе если:
   - заряды ядер этих атомов близки,
   - разница в теоретических значениях этих амплитуд не велика.

Общий вид команды MINIMIZE=GEDINT:

    MINIMIZE=GEDINT,N1-N2,N3-N4,.....,g1,g2,g3,....

Где:

- N1,N3 - номера категорий интенсивностей (например расстояния сопло-пластика),
- N2,N4 - номера кривых интенсивностей с соответствующих расстояний. Одним словом, N1,N2,N3....задают те кривые, на которых будет вестись минимизацию.
- g1,g2,g3 - номера групп параметров, которые будут варьироваться при выполнении данной команды MINIMIZE (не обязательно их указывать).

Количество кривых и номеров групп в команде MINIMIZE не ограничено.
Пример использования команды MINIMIZE:

    MINIMIZE=GEDINT,1-1

Что означает что минимизация идет по кривой sM(s), полученной из первой интенсивности первой категории. (1(категория)-1(номер)). Варьироваться будут все группы параметров, отличные от нуля.
Можно проводить минимизацию сразу по нескольким sM(s) (например по двум расстояниям сопло-пластинка сразу). В таком случае кривые после слова MINIMIZE указываются через пробел или запятую. Например:

    MINIMIZE=GEDINT,1-1,2-1

Что означает что минимизация одновременно пойдет по первой кривой первой категории и первой кривой второй категории.

    MINIMIZE=GEDINT,1-1,2-1,2,3,4,12,45

Минимизация функционала на кривых интенсивности 1-1 и 2-1, варьированием групп 2,3,4,12 и 45.
Следует обратить внимание на то что если Вы в команде MINIMIZE указываете номера групп варьируемых параметров, то если вы хотите варьировать и Scale-множители экспериментальных sM(s), - их номера групп тоже нужно указывать. Для этого лучше всего эти номера задать на этапе ввода интенсивностей (см. секцию ввода интенсивностей). Если же не задавать их, то UNEX самостоятельно присвоит им некоторые отличные от нуля числа, которые можно посмотреть в выходном файле. Однако такой вариант хуже: UNEX автоматически задает номера групп для Scale-множителей исходя из уже использующихся (для геометрии, амплитуд, ...) номеров групп. Так что если вы в какой-то момент добавите новые группы, например в амплитудах какой-либо молекулы, то UNEX может назначить параметрам Scale уже другие номера групп и Вам придется корректировать команду MINIMIZE.

## SEARCH

Поиск минимального значения функционала.
Для того чтобы искать минимальный R-фактор по параметру/параметрам нужно выполнить некоторые действия:

1. Получить экспериментальную кривую (кривые). Ее можно ввести в программу отдельно с помощью команды SMS или вычислить при помощи проведения линии фона (семейство команд BGL).
2. Запустить поиск минимального R-фактора с помощью команды SEARCH=метод,GEDINT...

Уже указывалось что результаты работы команды SEARCH особо важны при предположительном наличии нескольких решиний обратной задачи электронографии. В этом случае, перед началом оценки параметров по методу найменьших квадратов, важно знать хорошее начальное приближение, чтобы команда MINIMIZE гарантированно пришла в нужный нам минимум. Часто молекулы с несколькими волчками (торсионные степени свободы с относительно низкими барьерами) могут давать такие интенсивности рассеяния электронов, при интерпретации которых можно получить различные решения. На рис. 2 показан пример такой ситуации, когда на поверхности R-фактора в зависимсоти от двух торсионных углов обнаруживается два минимума.
Рис. 2: Поверхность значений R-фактора в зависимости от двух торсионных углов. Минимумы обозначены кружками.
Такую поверхность (или гиперповерхность) можно получить с помощью команды SEARCH=SCAN, однако если мы хотим чтобы команда SEARCH варьировала одновременно много параметров (например у нас много торсионных углов или сложная смесь молекул и т.д.), то тогда имеет смысл использовать метод Монте-Карло по команде SEARCH=RAND. Например если мы хотим "просканировать" пять параметров, каждый в ста точках, то общее количество точек будет 10 миллиардов. А если использовать "рандомизацию" SEARCH=RAND, то можно гораздо быстрее найти глобальный минимум.

# Определение геометрии молекул по вращательным постоянным

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

1. Ввести в поле BASE имя изучаемой молекулы. Обязательно должно быть поле у этой молекулы.
2. В поле изучаемой молекулы указать величины вращательных постоянных (если они есть), их ошибки и поправки к вращательным постоянным.
3. Если у этой молекулы есть изотопомеры с известными вращательными постоянными, то их тоже нужно ввести с помощью ключа isotopmols. Эти изотопомеры с точки зрения UNEX - тоже обычные молекулы, так что для каждой из них должно быть свое поле, где в данном случае важно указать вращательные постоянные, их ошибки и т.д.
4. Ввести геометрии всех молекул. Это делается как обычно с помощью z-матриц.
5. Выполнить команду MINIMIZE.

Замечания:

- Если вводятся вращательные постоянные для нулевого колебательного уровня и значения поправок не вводятся - то полученная таким образом геометрия (например по команде MINIMIZE) будет соответствовать r(0)-структуре. Если вводятся поправки на равновесные вращательные постоянные - то полученная геометрия будет равновесной. Т.е. результирующий тип геометрии напрямую зависит от того какие вводятся вращательные постоянные и поправки.
- При вводе геометрий изотопомеров важно указывать массы изотопов в z-марицах. Для основных природных изотопов UNEX сама выберет массу. Для других изотопов массы нужно вводить самостоятельно. Проверить массы атомов в молекулах можно выполнив команду PRINT=XYZ,mol_name При этом в лог будет выдана таблица декартовых координат и масс атомов.
- Номера групп параметров можно указывать только в z-матрице главной изучаемой молекулы, а в z-матрицах изотопомеров не обязательно.
- Если Вы желаете чтобы геометрия у всех изотопомеров была одинаковая - значения параметров в z-матрицах должны быть одинаковыми. Если есть желание выйти за пределы приближения Борна-Оппенгеймера - то вы можете задавать разные значения параметров разным изотопомерам. При варьировании групп, разности между параметрами изотопомеров будут сохранены.
- Обратите внимание на знак поправок к вращательным постоянным. В UNEX эти поправки определяются так: d=B(e)-B(v).

## MINIMIZE

Минимизация функционала, построенного на вращательных постоянных молекулы и ее изотопомеров осуществляется по команде:

    MINIMIZE=MOLMW

MOLMW означает функционал:
Где: i - номер изотопомера молекулы, j - номер вращательной постоянной i-го изотопомера (A, B или C), w - вес (равен обратному квадрату ошибки), символы "e" и "t" обозначают эксперимент и теорию соответственно.
Все особенности команды MINIMIZE описаны в общей части (основах).
Так же как и при минимизации электронографического функционала в команде MINIMIZE=MOLMW можно самому указывать номера групп параметров, которые будут варьироваться. В данном случае варьировать можно только геометрические параметры.

# Совместная ГЭ+МВ обратная задача

Комбинированная задача, которой отвечает гибридный функционал:
осуществляется по команде:

    MINIMIZE=GEDINT+MOLMW,1-1,2-1

Особенности такие же, как и в случае использования индивидуальных функционалов GEDINT и MOLMW.