Application Center - Maplesoft

App Preview:

Optimization and Lagrange multipliers

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

lagrange.mws

Calculus with Maple

         by Dr. Jack Wagner,  Lehman College, City University of New York

>                                                            Copyright 2001

Optimization, Lagrange Multipliers

The method of Lagrange is multipliers explained and illustrated geometrically. concatenation operator

      The problem of optimization, just as in the single variable case, involves maximizing or minimizing some function subject to a set of constraints.  For ease of visualization we are going to work in R^3   with two constraint equations, but the ideas and methods are identical for higher dimensions and larger constraint sets.  Formal proofs of the method of Lagrange multipliers abound.  What we present here is a geometric picture and geometric explanation of the procedure.

Ex. 4.1

        We suppose that we are required to optimize   F(x,y,z)+y+2*z-x^2      subject to two constraints: F is restricted to the sphere, x^2+y^2+z^2 = 25    and  also to the cylinder,   x^2+y^2 = 9 .
   
      
 First, examine a plot of the two constraints together with several level sets of F.  

>    restart: with(plots):with(linalg):

Warning, new definition for norm

Warning, new definition for trace

>    F:=y+2*z-x^2:

       Now we create plot structures for several level sets of F.  The multiple 5 in the expression 5*k is necessary to separate the sets by a sufficient amount to assure visibility.

>    P1:=k->implicitplot3d(F=5*k,x=-10..10,y=-10..10,z=-10..10,axes=boxed):

>    S1:=seq(P1(k),k=2..4):

>    C1:=x^2+y^2+z^2=25:

Now we define the constraints and create plot structures for their graphs.

>    P2:=implicitplot3d(C1,x=-10..10,y=-10..10,
z=-10..10,axes=boxed,color=red,style=wireframe):

>    C2:=x^2+y^2=9:

>    P3:=implicitplot3d(C2,x=-10..10,y=-10..10,
z=-10..10,axes=boxed,color=red,style=wireframe):

By inspection of the constraint equations they intersect where z = 4. This intersection is plotted as a spacecurve.  

>    S2:=spacecurve({[x,sqrt(9-x^2),4],[x,-sqrt(9-x^2),4],[x,sqrt(9-x^2),-4],[x,-sqrt(9-x^2),-4]},x=-3..3,numpoints=1000,color=black,thickness=2):

>    display(S1,P2,P3,S2);

[Maple Plot]

The surfaces of C1 and C2 intersect in two circles, S2.  At an optimal point for F, one of the level sets must be tangent to one of the circles in S2.  Looking ahead we reproduce the plot of the level set, F^(-1) (11)  , the maximum of F restricted to C1 and C2, together with the tangents to S2 and   F^(-1) (11) , and the gradients of F , C1 and C2 as well as the normal plane to S2 at the point of tangency.

>    F:=y+2*z-x^2:           #Function to be optimized.
C1:=x^2+y^2+z^2=25:      #Constraint
C2:=x^2+y^2=9:           #Constraint
Sa:=spacecurve({[x,sqrt(9-x^2),4],[x,-sqrt(9-x^2),4]},x=-3..3,numpoints=1000,color=black,thickness=2):  # Plot of intersection of C1 and C2
e:=[x,sqrt(9-x^2),4]:       # Positive part of intersection of C1 and C2
`de/dx`:=map(diff,e,x):
`de/dx`:=subs(x=0,`de/dx`):
T:=evalm([0,3,4]+s*eval(`de/dx`)): #Tangent to e.
Sb:=spacecurve(T,s=-4..4,color=red,thickness=3):
Sc:=implicitplot3d(F=11,x=-3..3,y=-3..5,z=0..6,color=blue,style=wireframe):
                             # Plot of level set of inverse image of 11.
r:=[x,y,z]:
gF:=grad(F,r):
gF:=subs(x=0,eval(gF)):
G1:=evalm([0,3,4]+s*gF):   #Gradient vector of F at P=[0,3,4]
Sd:=spacecurve(G1,s=-1.5..1.5,color=black,thickness=3):
gC1:=grad(lhs(C1),r):
gC1:=subs({x=0,y=3,z=4},eval(gC1)):
gC2:=grad(lhs(C2),r):
gC2:=subs({x=0,y=3,z=4},eval(gC2)):
G2:=evalm([0,3,4]+s*gC1):        #Gradient vector of C1 at P.
G3:=evalm([0,3,4]+s*gC2):       #Gradient vector of C2 at P.
G4:=evalm([0,3,4]+s*gC1+t*gC2):   #Plane determined by G3 and G2
Se:=spacecurve({G2,G3},s=-0.3..0.4,color=red,thickness=3):
Sf:=plot3d(G4,s=-.3..0.4,t=-.3..0.4,color=green,style=wireframe):
G5:=innerprod(evalm((r-[0,3,4])),gF)=0: #Plane orthogonal to G1
Sg:=implicitplot3d(G5,x=-2..2,y=-1..6,z=2..6,color=pink,style=wireframe):
display(Sa,Sb,Sc,Sd,Se,Sf,Sg,axes=framed,scaling=constrained);

[Maple Plot]

Referring to this illustration, we will identify each function, its graph and its plot by the same name.  

By the implicit function theorem, the level sets of F,   {   F^(-1) (k) } , yield surfaces. To satisfy the constraints, these are required to intersect the graphs of C1 and C2. It is clear that to satisfy both constraints simultaneously, the graphs of the level sets must  intersect their intersection, one of the two circles, S2.  In addition, it is  evident that the extreme values of F will occur at those points where a level set of F is tangent to the graph of  S2. This is analogous to the situation in linear programming where the optimized linear function must be tangent to the feasible set.  Call such a point of tangency P. At every point the tangent line to S2 must lie in the tangent plane to the graphs of both C1 and C2, and, therefore, in the intersection of those tangent planes.  In particular, this is true at P. Since, at P there is a graph of   {   F^(-1) (k) }  tangent to S2, for some value of k, say   k[0] ,   at this point the tangent line to S2 lies in the tangent plane to the graph of   { F^(-1) ( k[0] ) }   . At each point on C1 and C2, its gradient is orthogonal to its tangent plane.  Because the tangent line to S2 lies in the intersection of the tangent planes of C1 and C2, both gradients must be orthogonal to the tangent line to S2.  That is, at P on S2, the gradient vectors to C1 and C2 lie in a plane orthogonal to the tangent to S2, i.e. the normal plane to S2 at P.  Since the gradient vectors are independent any vector in that plane is a linear combination of the two gradients.  But the gradient of F, at P, is orthogonal to its tangent plane at that point, i.e. the plane tangent to the graph of   { F^(-1) ( k[0] ) }   . At P where F is just tangent to S2, the tangent plane of F must contain the tangent line to S2.  The gradient of F, therefore, is orthogonal to that tangent line and  so lies in the normal plane at P. The result of this is that the gradient of F is a linear combination of the gradients of the two constraint surface equations.      grad(F) = lambda[1]*grad(C[1])+lambda[2]*grad(C[2])    This yields three equations in five unknowns.  The lambdas are the Lagrange multipliers.  Together with the two constraint equations, we have five equations in five unknowns.  Implementing this in Maple is a lot simpler than explaining it.
 

>    r:=[x,y,z]:

>    gF:=grad(F,r);

gF := vector([-2*x, 1, 2])

>    gC:=evalm(lambda[1]*grad(lhs(C1),r)+lambda[2]*grad(lhs(C2),r));  

gC := vector([2*lambda[1]*x+2*lambda[2]*x, 2*lambda[1]*y+2*lambda[2]*y, 2*lambda[1]*z])

>    for k to 3 do eq.k:=gF[k]=gC[k] od;       

eq1 := -2*x = 2*lambda[1]*x+2*lambda[2]*x

eq2 := 1 = 2*lambda[1]*y+2*lambda[2]*y

eq3 := 2 = 2*lambda[1]*z

Note the use of the concatenation operator ".".

sol:=solve({eq1,eq2,eq3,C1,C2},{x,y,z,lambda[1],lambda[2]});

sol := {lambda[2] = -1/12, y = 3, z = 4, lambda[1] = 1/4, x = 0}, {lambda[2] = 5/12, y = 3, z = -4, lambda[1] = -1/4, x = 0}, {lambda[2] = -5/12, y = -3, z = 4, lambda[1] = 1/4, x = 0}, {lambda[2] = 1/...
sol := {lambda[2] = -1/12, y = 3, z = 4, lambda[1] = 1/4, x = 0}, {lambda[2] = 5/12, y = 3, z = -4, lambda[1] = -1/4, x = 0}, {lambda[2] = -5/12, y = -3, z = 4, lambda[1] = 1/4, x = 0}, {lambda[2] = 1/...
sol := {lambda[2] = -1/12, y = 3, z = 4, lambda[1] = 1/4, x = 0}, {lambda[2] = 5/12, y = 3, z = -4, lambda[1] = -1/4, x = 0}, {lambda[2] = -5/12, y = -3, z = 4, lambda[1] = 1/4, x = 0}, {lambda[2] = 1/...

Now we substitute each of these sets of values in turn into F.

>    V:=evalf(seq(subs(sol[k],F),k=1..6));

V := 11., -5., 5., -11., -1.250000000, -17.25000000

From whence it is clear that F is maximized at [0,3,4] and minimized at [`+/-`*sqrt(35)/2, -1/2, -4]

In order to understand this very important method more thoroughly we will dissect a second example.  

Ex. 4.2

Find the maximum and minimum values of the function: [Maple OLE 2.0 Object]

subject to the two constraints,   [Maple OLE 2.0 Object] and [Maple OLE 2.0 Object]

We will start by subjecting H to only the first constraint, c1.

restart: with(plots):with(linalg):

Warning, new definition for norm

Warning, new definition for trace

>    c1:=x^2+y^2-z=0: c2:=x^2+3*y^2=1:  H:=5*x^2+y^2+5*z^2-6*z:

>    r:=[x,y,z]:

>    `grad(H)`:=grad(H,r);

`grad(H)` := vector([10*x, 2*y, 10*z-6])

>    `grad(c1)`:=grad(lhs(c1),r);

`grad(c1)` := vector([2*x, 2*y, -1])

Set up the three equations containing   lambda[1] .

>    eq:=seq(`grad(H)`[k]=evalm(lambda[1]*`grad(c1)`[k]),k=1..3);

eq := 10*x = 2*lambda[1]*x, 2*y = 2*lambda[1]*y, 10*z-6 = -lambda[1]

Solve the system of equations consisting of the constraint equation and the three equations containing lambda[1] .

Note how the use of backward quotes on the left of the assignment operator permits the creation of more suggestive notation. The subscript selector must be placed outside the quotes.

>    sol_1:=solve({eq,c1},{x,y,z,lambda[1]});

sol_1 := {x = 0, y = 0, lambda[1] = 6, z = 0}, {x = RootOf(-1+10*_Z^2), y = 0, lambda[1] = 5, z = 1/10}, {x = 0, y = RootOf(-1+2*_Z^2), z = 1/2, lambda[1] = 1}

Now we compute values of H for each of these solutions.

>    V1:=evalf(seq(subs(sol_1[k],H),k=1..3));

V1 := 0, -.500000000e-1, -1.250000000

We note a maximum value for H of 0 and a minimum of -1.25.  Compare the graphs of    H^(-1) (0)   and   H^(-1) (-1.25)   with the graph of c1.

>    `P(H)(0)`:=implicitplot3d(H=0,x=-1..1,y=
-1..1,z=0..2,axes=framed,numpoints=2000,color=blue,style=wireframe):

>    `P(c1)`:=implicitplot3d(c1,x=-1..1,y=-1..1,z=0..2,axes=framed,color=turquoise):

>    display(`P(H)(0)`,`P(c1)`);

[Maple Plot]

You will be able to convince yourself,   by rotating the plot on the screen, that H is tangent to c1 at its lowermost extreme. Nonetheless, this graph does not represent a solution to the optimization problem. There are many points on c1 that yield larger values of H than 0. An obvious choice is [1,1,2].

>    subs({x=1,y=1,z=2},H);

14

This is just one of many instances in which a solution to an equation is not necessarily thge solution to the problem. Remember, we set up equations based on conditions the solution must meet. It's the "if" side of a proposition. If X is a solution then it will satisfy this equation. However, the "only if" side, "If X is a solution to the equation then it is a solution to the problem." may not be true.

Now the same thing for H= -1.25

>    `P(H)(-1.25)`:=implicitplot3d(H=-1.25,x=
-1..1,y=-1..1,z=0..2,axes=framed,numpoints=2000,color=blue,style=wireframe):

>    `P(c1)`:=implicitplot3d(c1,x=-1..1,y=-1..1,z=0..2,axes=framed,color=turquoise):

>    display(`P(H)(-1.25)`,`P(c1)`);

[Maple Plot]

Here the graph od H is tangent to the graph of c1 at each end of the ovoid graph of H. This is  a solution. There are no points on c1 that will make the value of H smaller than -1.25.

Now we turn our attention to optimizing H on c2 alone.

>    `grad(c2)`:=grad(lhs(c2),r);

>   

`grad(c2)` := vector([2*x, 6*y, 0])

Set up the equations containing lambda[2] .

>    eq2:=seq(`grad(H)`[k]=evalm(lambda[2]*`grad(c2)`[k]),k=1..3);

eq2 := 10*x = 2*lambda[2]*x, 2*y = 6*lambda[2]*y, 10*z-6 = 0

Solve the system consisting of the constraint and the equations containing   lambda[2] .

>    sol_2:=solve({eq2,c2},{x,y,z,lambda[2]});

sol_2 := {x = 0, z = 3/5, lambda[2] = 1/3, y = RootOf(-1+3*_Z^2)}, {y = 0, x = 1, z = 3/5, lambda[2] = 5}, {y = 0, z = 3/5, lambda[2] = 5, x = -1}

Compute the value of H at each point.

>    V2:=evalf(seq(subs(sol_2[k],H),k=1..3));

V2 := -1.466666667, 3.200000000, 3.200000000

Now we compare the plot of H=3.2 to the plot of c2.

>    `P(H)(3.2)`:=implicitplot3d(H=3.2,x=-1..1,
y=-1..1,z=0..2,axes=framed,numpoints=2000,color=blue):

>    `P(c2)`:=implicitplot3d(c2,x=-1..1,
y=-1..1,z=0..2,axes=framed,color=red,style=wireframe):

>    display(`P(H)(3.2)`,`P(c2)`);

[Maple Plot]

H is tangent to c2 at its outermost extreme. Again, this plot does not represent a solution to the optimization problem. It is certainly clear that many points on c2 will yield values for H larger than 3.2.  

>    P(H)(-1.466666667):=implicitplot3d(H=-1.466666667,x=-1..1,
y=-1..1,z=0..2,axes=boxed,numpoints=2000,color=blue):

>    display(P(H)(-1.466666667),`P(c2)`);

[Maple Plot]

Now the graph of H fits just inside the graph of c2.  This is a minimum. There is no level set   [Maple OLE 2.0 Object]

  | p<-1.466666667, which intersects the graph of c2.

Now we put it all together.  Since both constraints apply, we need to look at the intersection of c1 and c2.

>    E:=eliminate({c1,c2},z);

E := [{z = x^2+y^2}, {x^2+3*y^2-1}]

Here we have z in terms of x and y. Now we need y in terms of x.

>    solve(E[2,1],y);

1/3*sqrt(-3*x^2+3), -1/3*sqrt(-3*x^2+3)

This gives two parameterizations which, together, cover the line of intersection.  

>    L1:=[t,1/3*sqrt(-3*t^2+3),t^2+(1/3*sqrt(-3*t^2+3))^2]:

>    L2:=[t,-1/3*sqrt(-3*t^2+3),t^2+(1/3*sqrt(-3*t^2+3))^2]:

>    spacecurve({L2,L1},t=
-3..3,color=black,thickness=3,numpoints=1000,axes=framed);

[Maple Plot]

The space is an artifact of the plotting process.  We used 1,000 points to minimize this artifact.  

Now we will optimize H subject to both constraints.  Set up the three equations containing   lambda[1]   and   lambda[2] .

>    eq3:=seq(`grad(H)`[k]=evalm(lambda[1]*`grad(c1)`[k]+lambda[2]*`grad(c2)`[k]),k=1..3);

eq3 := 10*x = 2*lambda[1]*x+2*lambda[2]*x, 2*y = 2*lambda[1]*y+6*lambda[2]*y, 10*z-6 = -lambda[1]

Solve the system of five equations consisting of the two constraints and the three equations containing   lambda[1]   and   lambda[2] .

>    sol_3:=solve({eq3,c1,c2},{x,y,z,lambda[1],lambda[2]});

sol_3 := {x = 0, lambda[1] = 8/3, z = 1/3, y = RootOf(-1+3*_Z^2), lambda[2] = -5/9}, {y = 0, x = 1, lambda[2] = 9, lambda[1] = -4, z = 1}, {y = 0, lambda[2] = 9, lambda[1] = -4, z = 1, x = -1}, {y = 1/...
sol_3 := {x = 0, lambda[1] = 8/3, z = 1/3, y = RootOf(-1+3*_Z^2), lambda[2] = -5/9}, {y = 0, x = 1, lambda[2] = 9, lambda[1] = -4, z = 1}, {y = 0, lambda[2] = 9, lambda[1] = -4, z = 1, x = -1}, {y = 1/...

Here are the values of H corresponding to each of these solutions.

>    V3:=evalf(seq(subs(sol_3[k],H),k=1..4));

V3 := -1.111111111, 4., 4., -2.050000001

Since the largest value of H is 4 we plot   [Maple OLE 2.0 Object]  together with the curve of inteersection of the two constraints.

>    P(H)(4):=implicitplot3d(H=4,x=-2..2,y=-2..2,
z=-2..2,axes=framed,numpoints=2000,color=turquoise,style=wireframe):

>    Const:=spacecurve({L2,L1},t=
-3..3,color=black,thickness=3,numpoints=1000,axes=framed):

>    display(P(H)(4),Const);

[Maple Plot]

This is the plot of the maximum H.  

The sol_3 corresponding to -2.050000001, is imaginary. The smallest value of H corresponding to a real solution is, -1.111111111.  Plot   [Maple OLE 2.0 Object]

>    P(H)(-1.111111111):=implicitplot3d(H=-1.111111111,x=-1..1,
y=-1..1,z=0..1.5,axes=framed,numpoints=2000,color=turquoise,style=wireframe):

>    display(Const,P(H)(-1.111111111));

[Maple Plot]

>   

This is the plot of the minimum H.

      It is clear that restriction of H by two constraints is the same as restriction by a single constraint represented by the line of intersection of the two graphs representing the constraint equations.  

      The power and usefulness of Maple in presenting these plots is apparent.  The question as to whether a particular solution really represents a maximum or minimum is frequently answerable by inspection that can then be verified.

Practice

1.What are the radii of the largest and smallest  spheres, centered at the origin,  
tangent to the ellipse described by
[Maple OLE 2.0 Object]  ?

 2. Use the method of Lagrange multipliers to find the radius of the sphere, centered
at (a,b) and tangent to the line Ax+By+Cz=K, and thereby derive a formula for
the distance from a point to a line.

3. For each of the following functions, and sets of constraints:
      a) Plot the sets of constraints.
      b) Find the maximum and minimum of the functions subject to the indicated constraints.
      c) Plot the inverse image of the solutions together with the constraints.

                     1. [Maple OLE 2.0 Object]              10*z-5*x^2-3*y^2 = 20

                                                                   z+3*x^2+y^2 = 20

                     2. [Maple OLE 2.0 Object]                       (x-z)^2+(z-2)^2 = 1  

                                                                   (y-z)^2+(z-2)^2 = 1