adapted from http://en.wikipedia.org/wiki/Lagrange_multipliers.
clear all; close all
Suppose we seek to maximize the function subject to the constraint that . The function we seek to maximize is an unbounded plane, while the constraint is a unit circle. We plot these two functions here.
x = linspace(-1.5,1.5); [X,Y] = meshgrid(x,x); figure; hold all surf(X,Y,X+Y); shading interp;
coordinates, then convert those coordinates to x,y cartesian coordinates. then we plot f(x1,y1).
theta = linspace(0,2*pi); [x1,y1] = pol2cart(theta,1); plot3(x1,y1,x1+y1) view(38,8) % adjust view so it is easy to see
To find the maximum, we construct the following function: where , which is the constraint function. Since , we aren't really changing the original function, provided that the constraint is met!
function gamma = func(X) x = X(1); y = X(2); lambda = X(3); gamma = x + y + lambda*(x^2 + y^2 - 1); end
The minima/maxima of the augmented function are located where all of the partial derivatives of the augmented function are equal to zero, i.e. , , and . the process for solving this is usually to analytically evaluate the partial derivatives, and then solve the unconstrained resulting equations, which may be nonlinear.
Rather than perform the analytical differentiation, here we develop a way to numerically approximate the partial derivates.
function dLambda = dfunc(X)
function to compute the partial derivatives of the function defined above at the point X. We use a central finite difference approach to compute .
% initialize the partial derivative vector so it has the same shape % as the input X dLambda = nan(size(X)); h = 1e-3; % this is the step size used in the finite difference. for i=1:numel(X) dX=zeros(size(X)); dX(i) = h; dLambda(i) = (func(X+dX)-func(X-dX))/(2*h); end
The function we defined above (dfunc) will equal zero at a maximum or minimum. It turns out there are two solutions to this problem, but only one of them is the maximum value. Which solution you get depends on the initial guess provided to the solver. Here we have to use some judgement to identify the maximum.
X1 = fsolve(@dfunc,[1 1 0],optimset('display','off')) % let's check the value of the function at this solution fval1 = func(X1) sprintf('%f',fval1)
I don't know why the output in publish doesn't show up in the right place. sometimes that happens. The output from the cell above shows up at the end of the output for some reason.
Lagrange multipliers are a useful way to solve optimization problems with equality constraints. The finite difference approach used to approximate the partial derivatives is handy in the sense that we don't have to do the calculus to get the analytical derivatives. It is only an approximation to the partial derivatives though, and could be problematic for some problems. An adaptive method to more accurately estimate the derivatives would be helpful, but could be computationally expensive. The significant advantage of the numerical derivatives is that it works no matter how many variables there are without even changing the code!
There are other approaches to solving this kind of equation in Matlab, notably the use of fmincon.
ans = done
end % categories: optimization
X1 = 0.7071 0.7071 -0.7071 fval1 = 1.4142 ans = 1.414214