Добрый день!
Я занимаюсь научными исследованиями в области систем управления, и Matlab — мой основной рабочий инструмент. Одна из возможностей в MatLab — численная оптимизация. Оптимизировать (минимизировать) можно любую функцию, которая принимает на вход вектор варьируемых параметров и возвращает значение минимизируемого критерия. Естественно, в процессе оптимизации целевая функция вызывается множество раз и ее быстродействие существенно. В матлабе есть хорошие программные средства, которые часто позволяют существенно улучшить быстродействие, сохранив при этом читаемость и удобство сопровождения кода. Я приведу пример задачи, покажу на нём, что такое anonymous functions и nested functions, а потом покажу, как можно совместить эти два инструмента для заметного повышения быстродействия.
Постановка задачи и оптимизация
Я долго думал над примером для оптимизации, пытался выбрать что-то реалистичное – поиск оптимальной траектории для кинематически избыточного манипулятора, аппроксимация экспериментальных данных, идентификация параметров. Но все это как-то уводило от сути, так что решил остановится на непрактичном, но иллюстративном примере. Итак, заданы параметры a и b. Надо найти x в диапазоне [0;1], максимизирующий критерий:
где — некоторая функция, не зависящая явно от x, но нужная для вычисления критерия. Пусть это будет модуль суммы первого миллиона псевдослучайных чисел, полученных при установка seed в z. Напишем матлабовскую функцию, вычисляющую значение критерия:
function J = GetJ(a,b,x)
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));
J=-(phi_a*x^2+phi_b*sqrt(1-x^2));
end
Обратите внимание, что возвращаем значение со знаком минус, так как мы хотим найти максимум, а оптимизация ищет минимум.
Теперь для оптимизации при конкретных значениях a и b нам надо сделать анонимную функцию. Делается это так:
a=121233;
b=31235;
fhnd_GetJ = @(x) GetJ(a,b,x);
Теперь переменная fhnd_GetJ содержит хендл анонимной функции, которая принимает один параметр и позволяет вычислить GetJ(a,b, принятый параметр) для значений a и b, указанных при создании этой анонимной функции. Можно перейти непосредственно к минимизации. Скажем, мы хотим найти оптимальное значение x с точностью до 1 милионной:
opt=optimset('TolX',1e-6);
optimal_x=fminbnd(fhnd_GetJ,0,1,opt);
Функция fminbnd(fun,x_min,x_max) минимизирует скалярную функцию скалярного аргумента на интервале [x_min; x_max]. Здесь fun — хендл оптимизируемой функции. При выполнении оптимизации функция GetJ вызывается 12 раз, это можно увидеть, задав в опциях параметр ‘Display’ как ‘iter’ – покажет все итерации. На моем компе оптимизация занимает, в среднем, 10 мс.
Повышение быстродействия
Посмотрим, как это можно оптимизировать. Очевидно, что каждый раз при вызове функции GetJ мы совершенно зря считаем значения phi_a и phi_b, так как они при вариации x не меняются и зависят только от заданных значений a и b. Как тут можно сэкономить? Самый частый вариант, который мне встречается – вынести предварительные расчеты из целевой функции. Сделать вот такую функцию
function J = GetJ2(phi_a,phi_b,x)
J=-(phi_a*x^2+phi_b*sqrt(1-x^2));
end
и вот такой код:
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));
fhnd_GetJ2 = @(x) GetJ2(phi_a,phi_b,x);
optimal_x=fminbnd(fhnd_GetJ2,0,1,opt);
Конечно, это считается намного быстрее. Оптимизация проходит, в среднем, за 0.8 мс. Но платой за это является ухудшение кода. Разрывается целостность функционала – вычисление phi_a и phi_b самостоятельной ценности не имеет и в отрыве от вычисления J не нужно. Если где-то в другом месте опять потребуется считать J(a,b,x), уже не для оптимизации, а просто так, то придется вместо простого вызова GetJ за собой таскать еще и вычисление phi_a и phi_b или делать отдельно функции для оптимизации, отдельно для вычислений. Ну и просто не очень красиво.
Хорошо, что у матлаба есть способ обойти эту ситуацию. Для этого давайте создадим nested function внутри функции GetJ:
function J = GetJ3(a,b,x)
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));
J=nf_GetJ(x);
function out_val = nf_GetJ(x)
out_val=-(phi_a*x^2+phi_b*sqrt(1-x^2));
end
end
Nested function nf_GetJ видит все переменные в области видимости родительской функции и, в принципе, понятно, что и как она делает. Пока что мы никакого выигрыша по быстродействию не получили – все те же 10 мс на оптимизацию.
А теперь самое интересное – матлаб умеет работать с anonymous nested function. Вместо конкретного значения наша функция GetJ может вернуть хендл на собственную вложенную функцию. И при вызове функции по этому хендлу родительская функция выполняться не будет, но ее область видимости с посчитанными параметрами останется! Запишем:
function fhnd_J = GetJ4(a,b,x)
randn('seed',a);
phi_a=abs(sum(randn(1,1e6)));
randn('seed',b);
phi_b=abs(sum(randn(1,1e6)));
fhnd_J=@(x) nf_GetJ(x);
function out_val = nf_GetJ(x)
out_val=-(phi_a*x^2+phi_b*sqrt(1-x^2));
end
end
Теперь в возвращаемую переменную fhnd_J записывается хендл функции, позволяющей вычислить значение J для использованных параметров a и b, не вычисляя при этом переменные phi_a и phi_b, а используя значения, посчитанные при создании хендла. Далее наша оптимизация выглядит вот так:
fhnd_GetJ4=GetJ4(a,b,x);
optimal_x=fminbnd(fhnd_GetJ4,0,1,opt);
И выполняется такая оптимизация за 0,8 мс.
Результат
Мы получили такое же быстродействие, как при явном выносе вычислений в примере с GetJ2, но сохранили целостность функции и удобство ее использования.
В дальнейшем, если предполагается использовать функцию и для оптимизации, и просто для вычислений, то к ней можно добавить один необязательный параметр-флаг. Если он не задан, то функция возвращает значение. Если задан, то возвращается хендл на nested function.
Автор: Arastas
