Lagrange Multipliers in 2D

Set up function and solve for the constrained critical points

(( p * s + q * t) - (400000+195s+225t))/.{p->339-0.01s-.003t,q->399-0.004s-.01t}

Expand[%]

Clear[f,g1,g2,g3,ft,fs,
g1t,g1s,g2t,g2s,g3t,g3s,
s,t,L1,L2,L3,eqns,pts,gradf,gradg1,gradg2,gradg3]
f[t_,s_]:= -400000+144 s-0.01` (s^2)+174 t-0.007` s t-0.01` (t^2)
g1[t_,s_]:= s-5000
g2[t_,s_]:= t-8000
g3[t_,s_]:= s+t-10000

ft[t_,s_]=(∂_t) f[t,s]
fs[t_,s_]=(∂_s)f[t,s]
gradf[t_,s_]={ D[f[t,s],t], D[f[t,s],s]}

g1t[t_,s_]=(∂_t) g1[t,s]
g1s[t_,s_]=(∂_s)g1[t,s]
gradg1[t_,s_]={ D[g1[t,s],t], D[g1[t,s],s]}

g2t[t_,s_]=(∂_t) g2[t,s]
g2s[t_,s_]=(∂_s)g2[t,s]
gradg2[t_,s_]={ D[g2[t,s],t], D[g2[t,s],s]}

g3t[t_,s_]=(∂_t) g3[t,s]
g3s[t_,s_]=(∂_s)g3[t,s]
gradg3[t_,s_]={ D[g3[t,s],t], D[g3[t,s],s]}

eqns={
ft[t,s] == L1* g1t[t,s]+L2*g2t[t,s]+L3*g3t[t,s],
fs[t,s] == L1* g1s[t,s]+L2*g2s[t,s]+L3*g3s[t,s],
g1[t,s]==0,g2[t,s]==0,g3[t,s]==0}
pts=Solve[eqns,{t,s,L1,L2,L3}]

Clear[eqns2]
eqns2={
ft[t,s] == L3*g3t[t,s],
fs[t,s] == L3*g3s[t,s],
g3[t,s]==0}
pts=Solve[eqns2,{t,s,L3}]

Print["The method requires that the gradient of g is nonzero at a critical point"]
{gx[x,y], gy[x,y]}/.pts

f[x,y]/.pts
N[f[x,y]/.pts]

graphs

Clear[xlo,ylo,xhi,yhi,zlo,zhi,numCon,t1,t2];
xlo=ylo=-2; xhi=yhi=2;
numCon=25;

Clear[fPlot3D,fPlotGradient,fPlotLevel,
gPlot3Da,gPlot3Db,gPlot3DCurvea,gPlotGradient,gPlot3DCurveb,gPlotImplicit,
cpPlot2D,cpPlot3D];

fPlot3D=Plot3D[f[x,y],{x,xlo,xhi},{y,ylo,yhi},AxesLabel->{" x "," y "," z "}]

zlo=-1;zhi=20;

Solve[g[x,y]==0,y]

gPlot3Da=ParametricPlot3D[{x,-(Sqrt[1-(x^2)]),z},{x,xlo,xhi},{z,zlo,zhi}];
gPlot3Db=ParametricPlot3D[{x,(Sqrt[1-(x^2)]),z},{x,xlo,xhi},{z,zlo,zhi}];

Show[{fPlot3D,gPlot3Da,gPlot3Db},AxesLabel->{" x "," y "," z "}]

gPlot3DCurvea=ParametricPlot3D[{x,(Sqrt[1-(x^2)]),f[x,(Sqrt[1-(x^2)])]},{x,xlo,xhi},{z,zlo,zhi},PlotStyle->{Red,Thickness[.15]}];
gPlot3DCurveb=ParametricPlot3D[{x,-(Sqrt[1-(x^2)]),f[x,-(Sqrt[1-(x^2)])]},{x,xlo,xhi},{z,zlo,zhi},PlotStyle->{Red,Thickness[.15]}];
Show[{fPlot3D,gPlot3DCurvea,gPlot3DCurveb},AxesLabel->{" x "," y "," z "}]

Plot3D[f[x,y],{x,xlo,xhi},{y,ylo,yhi},AxesLabel->{" x "," y "," z "},RegionFunction->Function[{x,y,z},-0.1<g[x,y]<0.1]]

fPlotLevel=ContourPlot[f[x,y],{x,xlo,xhi},{y,ylo,yhi},Contours->numCon];
gPlotImplicit=ContourPlot[g[x,y]==0,{x,xlo,xhi},{y,ylo,yhi},Axes->True];
Show[{fPlotLevel,gPlotImplicit},Axes->True]

{Plot3D[f[x,y],{x,xlo,xhi},{y,ylo,yhi},AxesLabel->{" x "," y "," z "},ColorFunction->Hue],ContourPlot[f[x,y],{x,xlo,xhi},{y,ylo,yhi},Contours->numCon,ColorFunction->Hue]}

t1=Table[{x,y,f[x,y]}/.pts[[index]],{index,1,Length[pts]}]
cpPlot3D=ListPointPlot3D[t1,PlotStyle->PointSize[0.03]];
Show[{fPlot3D,cpPlot3D},AxesLabel->{" x "," y "," z "}]

Show[{Plot3D[f[x,y],{x,xlo,xhi},{y,ylo,yhi},AxesLabel->{" x "," y "," z "},RegionFunction->Function[{x,y,z},-0.1<g[x,y]<0.1]],cpPlot3D},AxesLabel->{" x "," y "," z "}]

t2=Table[Point[{x,y}/.pts[[index]]],{index,1,Length[pts]}];
cpPlot2D = Show[Graphics[{PointSize[0.02`],t2}],PlotRange->All,Axes->True];
Show[{{fPlotLevel,gPlotImplicit}},Graphics[{PointSize[0.02`],t2}],PlotRange->All,Axes->True]

CoolColor[ z_ ] :=RGBColor[z,1-z,1];
fPlotGradient=VectorPlot[gradf,{x,xlo,xhi},{y,ylo,yhi}, VectorStyle->Red]
gPlotGradient=VectorPlot[gradg,{x,xlo,xhi},{y,ylo,yhi},VectorStyle->Brown]

the critical points are on the constraint curve where the gradients of g and of f are parallel.

Show[{fPlotGradient,gPlotGradient,cpPlot2D,gPlotImplicit},PlotRange->All,Axes->True]

Show[{ContourPlot[f[x,y],{x,xlo,xhi},{y,ylo,yhi},Contours->numCon,ContourShading->False],fPlotGradient,gPlotGradient,cpPlot2D,gPlotImplicit},PlotRange->{{-1.25,1.25},{-1.25,1.25}},Axes->True]