Минимизация функции в Matlab

Question

Я хотел бы свести к минимуму (с использованием fmincon или аналогичный) следующую функцию:

function Difference= myfun3(wk,omega,lambda,Passetcovar,tau,PMat,i,Pi,Q) 
    wcalc=inv(lambda* Passetcovar)*inv(inv(tau * Passetcovar)+ PMat(i,:)'*inv(omega)*PMat(i,:))*(inv(tau * Passetcovar)*Pi+ PMat(i,:)'*inv(omega)*Q(i,:)); 
    Difference=sum((wk-wcalc).^2); 
end 

wk и wcalc являются < 8 х 1 двойные> векторы-столбцы, где WK известно и Wcalc дается выше уравнения.

Как минимизировать Difference варьируя Omega для Omega >0 с

• скалярной
лямбда
• Passetcovar- < 8 х 8 двойной>
• тау - Скалярное
• PMat- < 3 х 8 double>
• Omega- Scalar
• Q-3 х 1 двойной>
• < 8 Пи х 1 двойной>

Answer 1

Во-первых, sigma вектор-строка? Если нет, то f также является вектором. Вы пытаетесь многоцелевую оптимизацию? Тогда fminsearch не поможет.

Во-вторых, ознакомьтесь с документацией fminsearch перед ее использованием. Предполагается, что f является функцией, которая сопоставляет ваш входной вектор скаляру. Кроме того, ему нужна стартовая точка x0.

Таким образом, напишите функцию f, которая принимает omega и возвращает скалярное значение целевой функции. Кроме того, выясните допустимый x0 (то есть начальное значение для omega).

В-третьих, fminsearch не допускает ограничений. Вы можете взломать его, сделав f return Inf или что-то, когда omega <= 0. Я бы порекомендовал fmincon.

Ваша функция должна выглядеть следующим образом. Убедитесь, что все остальные переменные, такие как PMap,tau и т. Д. доступны по всему миру. В противном случае вам понадобится anonymous function, чтобы перейти к fminsearch.

obj = f(omega) 
wcalc=inv(lambda* sigma)*inv(inv(tau * sigma)+ PMap(i,:)'*inv(Omega)*PMap(i,:))*(inv(tau * sigma)*pi+ PMap(i,:)'*inv(Omega)*Q(i,:)); 
obj = sigma*(wk-wcalc).^2; 

Затем используйте fmincon. Предположим, что у вас есть начальное значение для omega.

fmincon(f,omega,[],[],[],[],0,Inf); 

The [] добавлены, так как мы хотим, чтобы ваши решения связаны снизу с помощью этой формы.

x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub) 

Действительно ли ваш f выглядит так?

obj = f(omega,PMap,sigma,.....) 

Где ...... представляет все другие переменные. Затем вы можете использовать anonymous functions следующим образом.

g = @(omega)f(omega,PMap,sigma,.....); 

Теперь вы можете использовать g в fmincon или fminsearch.

Answer 2

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

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

clear all 
close all 
clc 

f = @(x) (x+3).^2; 

x = linspace(-5,5,100); 
y = f(x); 
plot(x,y); 

ymin = fminsearch(f, 0); 
printf('Local min found at: %f\n', ymin); 

который отображает простой график и печатает:

 
Local min found at: -3.000000 

Обратите внимание, что вы должны укажите начальную точку для поиска. В этом случае я использовал 0. Когда вы только указываете fminsearch один параметр, он ожидает, что параметр будет структурой, которая не выглядит так, как вы используете.

От help fminsearch:

X = fminsearch(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a 
    structure with the function FUN in PROBLEM.objective, the start point 
    in PROBLEM.x0, the options structure in PROBLEM.options, and solver 
    name 'fminsearch' in PROBLEM.solver. The PROBLEM structure must have 
    all the fields. 

Скорее всего, вы хотите, чтобы это использование:

X = fminsearch(FUN,X0) starts at X0 and attempts to find a local minimizer 
    X of the function FUN. FUN is a function handle. FUN accepts input X and 
    returns a scalar function value F evaluated at X. X0 can be a scalar, vector 
    or matrix. 

Answer 3

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

Вот пример того, как это сделать (предполагая, что все необходимые переменные определены).

rangeMin = 0; 
rangeMax = 10; %Put an assumed upper bound for Omega here 
stepNumber = 10000; 
Omega = linspace(rangeMin, rangeMax, stepNumber); 
result = Inf(size(Omega)) 
for t = 1:length(Omega) 
    wcalc(t)=inv(lambda* sigma)*inv(inv(tau * sigma)+ PMap(i,:)'*inv(Omega(t))*PMap(i,:))*(inv(tau * sigma)*pi+ PMap(i,:)'*inv(Omega(t))*Q(i,:)); 
end 

wcalcMin = min(wcalc) 
OmegaMin = find(wcalc == wcalcMin) 

Если вы хотите увидеть, что функция на самом деле выглядит теперь вы можете сделать

plot(wcalc) 

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

Если вас устраивает поведение функции, вы можете повысить точность, установив rangeMin и rangeMax, чтобы оценить только интересующую вас область. Если вас не устраивает поведение, вы можете попробовать, помогает ли увеличение stepNumber.

Answer 4

Эта цель НЕ является ужасно сложной. Но вы написали это, чтобы выглядеть сложным.

Прежде всего, омега - это скаляр! Зачем писать inv (омега)? Разделение на omega - лучшая идея, так как для обратной функции не потребуются накладные расходы.

Далее, tau - известный постоянный скаляр, как и Passetcovar. Зачем вычислять обратную матрицу (tau * Passetcovar) каждый раз, когда вызывается функция? Не только это, но вы вычисляете ту же самую обратную матрицу три раза в одной строке. Учитесь прекомпомировать эти вещи. Это сэкономит вам много времени и много головных болей.

В любом случае, у вас есть навязчивая идея с inv. Он называется 6 раз в одной строке, и большинство из этих вызовов являются излишними.

Позволяет переписать эту единственную строчку. Прежде всего, прекомпуте inv (Passetcovar), передавая все это в вашу цель, поэтому вам нужно делать это ТОЛЬКО один раз.

Примечания основной идентичность:

inv(k*A) = inv(A)/k 

, которая справедлива для любых ненулевых скалярных к.

IP = inv(Passetcovar); 

Опять же, не повторить вычисление и (Passetcovar) каждый раз, когда целевая функция называется. Вместо этого, вычислите его ONCE, прежде чем вы начнете оптимизацию.

Таким образом, вычисление становится немного проще читать:

wcalc = IP./lambda*inv(IP./tau + PMat(i,:)'*PMat(i,:)/omega)*(IP*Pi./tau + PMat(i,:)'*Q(i,:)/omega); 

Edit:

Наконец мы узнаем, что Pi является вектором. Я предполагаю, что это должен быть вектор 8x1, чтобы совпадение массива соответствовало.

Мы можем сэкономить еще несколько разделов и умножить, разложив некоторые константы и вставив паренны в удобное место. Заметим, что, вычислив вектор матрицы *, умножьте FIRST, а затем умножьте на IP, мы преобразуем 8x8 * 8x8 умножить на умножение 8x8 X 8x1. Для небольших массивов, подобных этому, разница не огромна, но стоит вспомнить идею.

wcalc = IP*(inv(IP + tau*PMat(i,:)'*PMat(i,:)/omega)*(IP*Pi + PMat(i,:)'*Q(i,:)/omega*tau))/(lambda*tau^2); 

Цель состоит в том, чтобы минимизировать сумму квадратов между Wcalc и WK, при условии положительной омеги. Теперь это скаляр.

Я бы предложил сначала загладить функцию, чтобы узнать что-то о ее форме и посмотреть, что может быть хорошим стартовым значением для омеги. Таким образом, обернув дескриптор функции вокруг myfun, ezplot сделает сюжет красиво, здесь для диапазона [0,100] для омеги. Выберите свою собственную верхнюю границу для омеги, если это необоснованно.

ezplot(@(omega) myfun3(wk,omega,lambda,Passetcovar,tau,PMat,i,Pi,Q),[0,100]) 

Таким образом, тривиальное решение заключается в использовании fminbnd, обеспечивая некоторое разумное, но достаточно большое значение для верхней грани. Хорошая вещь о fminbnd - это не нужно начинать значения. Вам нужно будет выбрать разумную верхнюю границу для омеги. Дело в том, что используйте инструмент, предназначенный для минимизации скалярной функции. Общие оптимизаторы, такие как fmincon, не нужны и требуют начального значения.

finalomega = fminbnd(@(omega) myfun3(wk,omega,lambda,Passetcovar,tau,PMat,i,Pi,Q),[0,100]) 

Также вы можете использовать fminsearchbnd, найденные при обмене файлами. Он может свести к минимуму функцию, подчиненную только ограничениям, связанным с нижними границами, но для fminsearchbnd потребуется начальное значение для omega.