MATH 448/548 - Example 2.2. The Color TV Problem with Constraints (p.33-44)

The Lagrange Multiplier method is used here to solve a constraint optimization problem.

Contents

Introduce the symbolic (independent) variables

clear all; close all;
syms x1 x2 lambda

Define and plot the level sets of the objective function and the constaint.

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')

Apply the Lagrange multiplier method

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

Plot the optimal level set and optimal solution

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')

Check the values of y on the other boundaries of the feasible set

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');

Sensitivity analysis to paramenter *a*

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
 
 

Sensitivity to constraint parameter *c* (Shadow prices)

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