The Lagrange Multiplier method is used here to solve a constraint optimization problem.
clear all; close all; syms x1 x2 lambda
Note the use of the Matlab command ezcontourf. An alternate command for contour plots is ezcontour. Also note how one can add a colorbar and a customized title to your plot.
y = (339-.01*x1-.003*x2)*x1+(399-.004*x1-.01*x2)*x2-(400000+195*x1+225*x2); g = x1 + x2 - 10000; ezcontourf(y,[0 10000 0 10000]); hold on ezplot(g,[0 10000 0 10000]); colorbar('EastOutside') title('A few level sets of y and the constaint x_1+x_2 = 10000')
dydx1 = diff(y,x1);
dydx2 = diff(y,x2);
dgdx1 = diff(g,x1);
dgdx2 = diff(g,x2);
[lambdahat, x1hat, x2hat] = solve(dydx1-lambda*dgdx1, dydx2-lambda*dgdx2, g);
format bank
lambdahat=double(lambdahat)
x1hat = double(x1hat), x2hat = double(x2hat)
yhat=subs(y,[x1,x2],[x1hat,x2hat])
lambdahat = 24.00 x1hat = 3846.15 x2hat = 6153.85 yhat = 532307.69
ezcontourf(y,[0 10000 0 10000]); hold on; ezplot(y-yhat,[0 10000 0 10000]); hold on; ezplot(g,[0 10000 0 10000]); hold on; plot(x1hat,x2hat,'ko'); colorbar('EastOutside') title('Optimal level set y = y_{max} and constaint x_1+x_2 = 10000')
figure, ezcontourf(y,[0 5000 0 8000]); hold on; ezplot(y-yhat,[0 5000 0 8000]); hold on; ezplot(g,[0 5000 0 8000]); hold on; plot(x1hat,x2hat,'ko'); colorbar('EastOutside'); title('Optimal level set y = y_{max} and constaints x_1+x_2 \leq 10000, x_1 \leq 5000, x_2\leq 8000');
syms a y = (339-a*x1-.003*x2)*x1+(399-.004*x1-.01*x2)*x2-(400000+195*x1+225*x2); g = x1+x2-10000; dydx1 = diff(y,x1); dydx2 = diff(y,x2); dgdx1 = diff(g,x1); dgdx2 = diff(g,x2); [lambdahata, x1hata, x2hata] = solve(dydx1-lambda*dgdx1, dydx2-lambda*dgdx2, g) Sx1a=diff(x1hata,a)*a/x1hata; Sx1a=subs(Sx1a,a,0.01) Sx2a=diff(x2hata,a)*a/x2hata; Sx2a=subs(Sx2a,a,0.01) %see page 46 Sya=diff(y,a)*a/yhat; Sya=subs(Sya,a,0.01) %see page 49
lambdahata = -52*(500*a-11)/(3+1000*a) x1hata = 50000/(3+1000*a) x2hata = 20000*(-1+500*a)/(3+1000*a) Sx1a = -0.77 Sx2a = 0.48 Sya = -1073741824/57156103246769225*x1^2
syms c
y = (339-0.01*x1-.003*x2)*x1+(399-.004*x1-.01*x2)*x2-(400000+195*x1+225*x2);
g = x1+x2-c;
dydx1 = diff(y,x1);
dydx2 = diff(y,x2);
dgdx1 = diff(g,x1);
dgdx2 = diff(g,x2);
[lambdahata, x1hatc, x2hatc] = solve(dydx1-lambda*dgdx1, dydx2-lambda*dgdx2, g)
yhatc=subs(y,[x1,x2],[x1hatc,x2hatc])
dydc=diff(yhatc,c)
subs(dydc,c,10000)
lambdahata = 159-27/2000*c x1hatc = -15000/13+1/2*c x2hatc = 15000/13+1/2*c yhatc = (4512/13-13/2000*c)*(-15000/13+1/2*c)+(5097/13-7/1000*c)*(15000/13+1/2*c)-5650000/13-210*c dydc = 159-27/2000*c ans = 24.00