# Dowling13B.wxmx
Optimization with Inequality Constraints

load(draw)$  set_draw_defaults(line_width=2, draw_realpart=false)$ 
fpprintprec:5$   ratprint:false$     kill(all)$

load ("Econ1.mac");

## Preface

Dowling13B.wxmx uses Maxima to work some of the problems in  the last section of
Ch. 13 of Introduction to Mathematical Economics (3rd ed), by Edward T. Dowling, 
(Schaum's Outline Series),  McGraw-Hill, 2012. This text is a bargain, with many 
complete problems worked out in detail. You should compare Dowling's solutions, 
worked out "by hand", with what we do using  Maxima here.
 
Typical problems considered here are to maximize utility or profit while satisfying one or 
two constraints, which have the form g1(x,y) >= 0, g2(x,y) >= 0. 

Also considered is the problem of minimizing cost while satisfying one or two constraints, 
which have the form g1(x,y) <= 0, g2(x,y) <= 0.

Three methods are used for each maximization or minimization problem considered.

1). Purely graphical using a plot with lines and/or curves which represent the 
   objective function f(x,y) = constant, and the lines and/or curves which represent
   the constraints.

2). A grid search (which ignores points which do not satisfy the constraints) of 
    some (x,y) region, and

3). The KKT algorithm.

Solutions are checked by forming the Langrangian function and calculating
its first derivatives.

Two code files, Econ1.mac and KKT.mac are available in the same section (of Economic 
Analysis with Maxima), which define Maxima functions used in this worksheet.

Use load ("Econ1.mac") and load ("KKT.mac").

Edwin L. (Ted) Woollett
http://web.csulb.edu/~woollett/
July, 2021

## KKT(f, min/max, g ), KKT(f, min/max, g1, g2)

Maxima code KKT (Karush-Kuhn-Tucker). Use load ("KKT.mac") to load our 
definitions of the Maxima function KKT.

SYNTAX:

KKT (f(x,y), max,  g(x,y) ) with g(x,y) >= 0

KKT (f(x,y), max, g1(x,y), g2(x,y) ) with g1(x,y) >=0 and g2(x,y) >=0
------------------------
KKT (f(x,y),  min, g(x,y) ) with g(x,y) <= 0

KKT (f(x,y),  min, g1(x,y), g2(x,y) ) with g1(x,y) <=0 and g2(x,y) <=0
-------------------------

load ("KKT.mac");

## KKT Method Summary

References: 
1.Alpha C. Chiang and Kevin Wainwright, Fundamental Methods of
Mathematical Economics, 4th ed, McGraw-Hill, 2005 (Sec. 13.1).

2. https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions

For one constraint, say KKT(f, max, g), our code uses the first derivatives of the
Lagrangian function L = f + lam*g. Our code is limited to two variables which must
be called (x,y). Sticking to that limitation we then require that:

a solution to the maximization problem with one constraint satisfy
1.  ∂L/∂x <= 0, x >= 0, x*∂L/∂x = 0
2.  ∂L/∂y <= 0, y >= 0, y*∂L/∂y = 0
3.  ∂L/∂lam >= 0, lam >= 0, lam*∂L/∂lam = 0

a solution to the minimization problem with one constraint satisfy
1.  ∂L/∂x >= 0, x >= 0, x*∂L/∂x = 0
2.  ∂L/∂y >= 0, y >= 0, y*∂L/∂y = 0
3.  ∂L/∂lam <= 0, lam >= 0, lam*∂L/∂lam = 0.

Note that in each case above, ∂L/∂lam = g(x,y), so the conditions on ∂L/∂lam are the
same as the conditions on g.

For two constraints, say KKT(f, max, g1, g2), our code uses the first derivatives of the
Lagrangian function L = f + lam1*g1 + lam2*g2. Our code is limited to two variables 
which must be called (x,y). Sticking to that limitation we then require that:

a solution to the maximization problem with two constraints satisfy
1.  ∂L/∂x <= 0, x >= 0, x*∂L/∂x = 0
2.  ∂L/∂y <= 0, y >= 0, y*∂L/∂y = 0
3.  ∂L/∂lam1 >= 0, lam1 >= 0, lam1*∂L/∂lam1 = 0
4.  ∂L/∂lam2 >= 0, lam2 >= 0, lam2*∂L/∂lam2 = 0

a solution to the minimization problem with two constraints satisfy
1.  ∂L/∂x >= 0, x >= 0, x*∂L/∂x = 0
2.  ∂L/∂y >= 0, y >= 0, y*∂L/∂y = 0
3.  ∂L/∂lam1 <= 0, lam1 >= 0, lam1*∂L/∂lam1 = 0.
4.  ∂L/∂lam2 <= 0, lam2 >= 0, lam2*∂L/∂lam2 = 0.

Note that in each case above, ∂L/∂lam1 = g1(x,y), so the conditions on ∂L/∂lam1 are the
same as the conditions on g1 (and likewise for g2(x,y) and ∂L/∂lam2).

The extension to more than two variables and more than two constraints is then
straightforward.

## Dowling Prob. 13.35, Maximize Profit, 1 constraint

Maximize profits π = 64 x - 2 x^2 + 96 y - 4 y^2 - 13 subject to the production 
constraint x + y <= 20. 

If we define g(x,y) = 20 - x - y >= 0, then 20 >= x + y  ->   x + y  <=  20.

Let Pr stand for the profit.

Pr : 64*x - 2*x^2 + 96*y - 4*y^2 - 13;

### Graphical Search

Plot the "feasible set" of solution points (x,y).

wxdraw2d (yrange = [0,25], xrange = [0,25], key = " feasible set ", 
    filled_func = true, fill_color = light_gray, explicit (20 - x,x,0,30))$

Look at the values of profit Pr at the corners of the feasible set.

[at (Pr, [x = 0, y = 20]), at (Pr, [x = 20, y = 0]) ];

Plot some Pr = constant curves using the implicit function inside wxdraw2d.

wxdraw2d (yrange = [0,30], xrange = [0,30], line_width = 2, proportional_axes = 'xy,
    xlabel = "x", ylabel = "y",  key = " feasible set ", 
    filled_func = true, fill_color = light_gray, explicit (20 - x,x,0,30),
    filled_func = false, 
    color = dark_turquoise, key = "Pr = 800",  implicit (Pr = 800, x,3,18,y,4,18),
    color = red, key = "Pr = 900",  implicit (Pr = 900, x,3,15,y,4,18),
    color = blue, key = "Pr = 1010",  implicit (Pr = 1010, x,3,15,y,4,15))$

The plot implies the optimum (maximum) point (in the feasible set) will be found on the
boundary defined by x + y <= 20.

### KKT Search

Try using our function KKT

Pr;

g : 20  - x - y;

KKT (Pr,max, g );

KKT has found the optimum solution  with lam > 0 and hence g = 0, the constraint is 
binding; the maximum occurs on the boundary (of the feasible set) defined by g >=0, 
more specifically, the boundary defined by x + y <= 20.

Put the maximum value location on our figure using the function points inside wxdraw2d.

wxdraw2d (yrange = [0,30], xrange = [0,30], proportional_axes = 'xy,
     line_width = 2,  filled_func = true, fill_color = light_gray, explicit (20 - x,x,0,30),
    filled_func = false,
    color = black, key = "y = 9.33", explicit (28/3, x, 0, 32/3),
    color = purple, key = "x = 10.67", parametric (32/3, t, t, 0, 28/3),
    color = dark_turquoise, key = "Pr = 800",  implicit (Pr = 800, x,3,15,y,4,15),
    color = red, key = "Pr = 900",  implicit (Pr = 900, x,3,15,y,4,15),
    color = blue, key = "Pr = 1000",  implicit (Pr = 1000, x,3,15,y,4,15),
    color = black,
    point_type = filled_circle, key = "Pr = 989.67", points ([ [32/3, 28/3] ])  )$

To check the solution found, form the Lagrangian function L, and find the relevant
derivatives.

L : Pr + lam*g;
[Lx : diff(L,x), Ly : diff(L,y), Llam : diff (L, lam) ];

Then find their values at the solution found by KKT.

subst ([x=32/3,y=28/3,lam=64/3], [Lx, Ly, Llam,g]);

Since the optimum x, y, and λ values are all greater than zero, part of the KKT conditions 
for an optimum are that the corresponding derivatives of the Lagrangian function should 
be zero. So we can find the same solution using solve:

solve ([Lx,Ly,Llam], [x, y, lam] );

### gridSearch (func, xR, yR, maxmin, constraint)

In our software file Econ1.mac is a straightforward grid search function
we call gridSearch, with the following four acceptable invocations: (You could use
f(u,v), g(u,v), for example instead of f(x,y), g(x,y), or any other variable names.
The first symbol (x) in the second argument (list) is always taken as one of the two
variable names, and the first symbol (y) in the third argument (list) is always taken
as the second variable name.

gridSearch (f(x,y), [x, x1, x2, dx], [y, y1, y2, dy], max, g(x,y))
     returns [fmax, [xmax, ymax] ] , and only looks at points for which g(x,y) >= 0.
     
gridSearch (f(x,y), [x, x1, x2, dx], [y, y1, y2, dy], min, g(x,y))
     returns [fmin, [xmin, ymin] ] , and only looks at points for which g(x,y) <= 0.

gridSearch (f(x,y), [x, x1, x2, dx], [y, y1, y2, dy], max, g1(x,y), g2(x,y) )
     returns [fmax, [xmax, ymax] ] , and only looks at points for which both
             g1(x,y) >= 0  and g2(x,y) >= 0.
             
gridSearch (f(x,y), [x, x1, x2, dx], [y, y1, y2, dy], min, g1(x,y), g2(x,y) )
     returns [fmin, [xmin, ymin] ] , and only looks at points for which 
        both g1(x,y) <= 0  and g2 (x,y) <= 0.

gridSearch (Pr, [x,10,11,0.01], [y,9,10,0.01], max, g);

which is the same answer as we have been finding in other ways.

## Dowling Prob. 13.36, Continuation Problem

Maximize the production function in Prob. 13.35 with the new production constraint
  x + y <= 36, enlarging the feasible set. Let g(x,y) = 36 - x - y >= 0. Then 36 >= x + y.

Let's use the same wxdraw2d code we used with Prob. 13.35, and change the boundary
to x+y = 36. Add some curves of constant Pr using implicit. This takes some experimentation
with the values of the constant and the range of x and y used inside the implicit function.

wxdraw2d (yrange = [0,40], xrange = [0,40], xlabel = "x", ylabel  = "y",
     proportional_axes = 'xy,  line_width = 2, key = " feasible set ", 
    filled_func = true, fill_color = light_gray, explicit (36 - x,x,0,40),
    filled_func = false, 
    color = dark_turquoise, key = "Pr = 1000",  implicit (Pr = 1000, x,3,30,y,4,25),
    color = red, key = "Pr = 1050",  implicit (Pr = 1050, x,3,30,y,4,20),
    color = blue, key = "Pr = 1070",  implicit (Pr = 1070, x,3,30,y,4,20))$

If you create a separate gnuplot window by using draw2d instead of wxdraw2d:

draw2d (yrange = [0,40], xrange = [0,40], xlabel = "x", ylabel  = "y",
     proportional_axes = 'xy,  line_width = 2, key = " feasible set ", 
    filled_func = true, fill_color = light_gray, explicit (36 - x,x,0,40),
    filled_func = false, 
    color = dark_turquoise, key = "Pr = 1000",  implicit (Pr = 1000, x,3,30,y,4,25),
    color = red, key = "Pr = 1050",  implicit (Pr = 1050, x,3,30,y,4,20),
    color = blue, key = "Pr = 1070",  implicit (Pr = 1070, x,3,30,y,4,20));

you can expand the plot to full screen and then watch the point coordinates change
as you move the cursor around. Placing the cursor at the center of the blue oval, 
you have almost found by inspection that x ~ 16, y ~ 12. You have to close the
separate gnuplot window to continue your work in wxMaxima.

It is clear from this plot that the optimum (maximum) point will be "interior" to the feasible
set and the constraint is "not binding".

After defining a new expression for g >= 0, we can call gridSearch1 for an approximate
solution.

g : 36 - x - y;

gridSearch (Pr, [x, 15,17,0.1], [y, 11,13,0.1], max, g);

Next let's try out KKT, 

KKT (Pr, max, g);

wxdraw2d (yrange = [0,40], xrange = [0,40], xlabel = "x", ylabel  = "y",
    proportional_axes = 'xy, line_width = 2, key = " feasible set ", 
    filled_func = true, fill_color = light_gray, explicit (36 - x,x,0,40),
    filled_func = false, 
    color = dark_turquoise, key = "Pr = 1000",  implicit (Pr = 1000, x,3,30,y,4,25),
    color = red, key = "Pr = 1050",  implicit (Pr = 1050, x,3,30,y,4,20),
    color = blue, key = "Pr = 1070",  implicit (Pr = 1070, x,3,30,y,4,20),
    color = black, 
    point_type = filled_circle, key = "Pr = 1075", points ([ [16, 12] ]),
    key = "y = 12",   explicit (12, x, 0, 16), 
    color = purple, key = "x = 16", parametric (16, t, t, 0, 12) )$

Check solution using first derivatives of the Lagrangian function.

L : Pr + lam*g;
[Lx : diff(L,x),Ly : diff(L,y), Llam : diff (L, lam)];
subst ([x = 16,y = 12, lam = 0], [Pr,x + y,Lx, Ly, Llam,g]);

With x + y  = 28 < 36, the production constraint is met with a profit = 1075 and the
solution is optimal since none of the Kuhn-Tucker conditions are violated, and no other
possible scenario satisfies all the Kuhn-Tucker conditions. (The program KKT
systematically explores all possible scenarios.)

The solution returned by solve, IF we assume Lx = 0, Ly = 0, Llam = 0, is:

soln : solve ([Lx,Ly,Llam], [x,y,lam])[1];

and is SUB-OPTIMAL since λ =  - 64/3 < 0

subst (soln,[Pr,x + y,Lx, Ly, Llam,g]);

float(%);

This sub-optimal solution requires a greater production (x + y = 36) with less profit (989.67)
than the optimal solution.

## Dowling Prob. 13.37, Minimize Cost, one constraint

Dowling, PROB 13.37:

Minimize the cost c = 5 x^2 - 80 x + y^2 - 32 y subject to x + y >= 30.

We then need g(x,y) = 30 - x - y <= 0 so that 30 <= x + y. 
Note carefully: when calling KKT(f, min, g ),  we need g(x,y)  <= 0.

wxdraw2d (yrange = [0,60], xrange = [0,60],  
    key = "need x + y > 30",explicit (30 - x,x,0,30))$

Add some curves c = const to this plot.

c : 5*x^2 - 80*x + y^2 - 32*y;

wxdraw2d (yrange = [0,60], xrange = [0,60],  
    key = "need x + y >= 30",explicit (30 - x,x,0,30),
    color = dark_turquoise, key = "c = 100",  implicit (c = 100, x,0,40,y,0,40),
color = red, key = "c = 10",  implicit (c = 10, x,0,40,y,0,40),
color = brown, key = "c = -100",  implicit (c = -100, x,0,40,y,0,40),
    color = blue, key = "c = -400",  implicit (c = -400, x,0,40,y,0,40),
    color = magenta, key = "c = -500",  implicit (c = -500, x,0,40,y,0,40))$
    

Remember we are seeking a solution for which the cost(x,y) is minimum subject to
x + y >= 30,
From the plot we can find a minimum for which x + y = 30, and c < -500. In principle, we
could lock-in on this minimum by adding a few more cost curves. 

If you remove the wx from wxdraw for the above figure and interactively read the coordinates
as you move the cursor in the separate gnuplot window (see above), you can see that
roughly, x ~ 8.8, y ~ 21.6.

g : 30 - x - y;
gridSearch (c, [x, 8, 10, 0.1],[y, 20, 22, 0.1], min, g );

Next try KKT.

KKT( c, min, g );

wxdraw2d (yrange = [0,60], xrange = [0,60],  proportional_axes = 'xy,
    key = "need x + y >= 30",explicit (30 - x,x,0,30),
    color = dark_turquoise, key = "c = 100",  implicit (c = 100, x,0,40,y,0,40),
color = red, key = "c = 10",  implicit (c = 10, x,0,40,y,0,40),
color = brown, key = "c = -100",  implicit (c = -100, x,0,40,y,0,40),
    color = blue, key = "c = -400",  implicit (c = -400, x,0,40,y,0,40),
    color = magenta, key = "c = -500",  implicit (c = -500, x,0,40,y,0,40),
    color = black, point_type = filled_circle, 
    key = " (9,21)",   points( [ [9, 21] ] ))$

Solution analysis:

L : c + lam*g;
[Lx, Ly, Llam] : [diff(L,x), diff(L,y), diff(L,lam)];
subst ( [x = 9, y = 21,lam = 10], [c, x + y, Lx, Ly, g] );

We see this solution satisfies g <= 0, since g = 0.

## Dowling Prob. 13.38, Continuation Problem

Dowling Prob 13.38
Minimize the same cost function c = 5 x^2 - 80 x + y^2 - 32 y as in Prob 13.37 
subject to the new constraint x + y >= 20, enlarging the feasible set.
Now:   g(x,y) = 20 - x - y <= 0  =>  20 <= x + y.

wxdraw2d (yrange = [0,25], xrange = [0,25],  xlabel = "x",ylabel  = "y",
    key = "need x + y >= 20",explicit (20 - x,x,0,20),       
    color = magenta, key = "c = -500",  implicit (c = -500, x,0,40,y,0,40),
    color = dark_turquoise, key = "c = -550",  implicit (c = -550, x,0,30,y,0,30),
    color = blue, key = "c = -570",  implicit (c = -570, x,0,40,y,0,40))$

We see from the plot that we can find a minimum < -570 for which x + y > 20.
Using the draw2d version in a separate gnuplot window and watching the coordinates
of the cursor, we see roughly x ~ 8, y ~ 16.

Try using gridSearch (..., min, g) with the new g = 20 - x - y <= 0

g : 20 - x - y;
gridSearch (c, [x, 7, 9, 0.1], [y, 15,17, 0.1], min, g);

Next try KKT:

KKT( c, min, g );

With lam = 0, the constraint is nonbinding. The firm can minimize costs (using  the
optimal solution: x = 8 and y = 16) while exceeding its production quota of 20 (by 4 units).

Details of solution:

L : c + lam*g;
[Lx, Ly, Llam] : [diff(L,x), diff(L,y), diff(L,lam)];
subst ( [x = 8, y = 16,lam = 0], [c, x + y, Lx, Ly, g] );

We see that x + y = 24 > 20.

## Chris Jones KKT One Constraint Example

Chris Jones, 
Math. Dept, Univ. of Leeds
http://www1.maths.leeds.ac.uk/~cajones/pubs.html

Math 2640, Introduction to Optimization
http://www1.maths.leeds.ac.uk/~cajones/math2640/MATH2640.html

Sample KKT problem, pdf page 4:
http://www1.maths.leeds.ac.uk/~cajones/math2640/notes4.pdf

Maximize  f = x - y^2, s.t. x^2 + y^2 <= 4.  i.e., s.t. 0 <= 4 - x^2 - y^2.
So define g(x,y) = 4 - x^2 - y^2 >= 0.

f : x - y^2;
g : 4 - x^2 - y^2;

Make a plot of the feasible set of solution points (x,y) and 
curve of f(x,y) = constant.

wxdraw2d (yrange = [0, 3], xrange = [0,3], proportional_axes = 'xy,
     key = " need x^2 + y^2 <= 4 ", filled_func = true, fill_color = light_gray, 
     explicit (sqrt (4 - x^2), x, 0, 2), filled_func = false, 
    point_size  = 1.5, point_type = filled_circle,
    color = blue, key = "f = -4", points ([ [0, 2] ]),
    color = red, key = "f = 0",  implicit (f = 0, x,0,2,y,0,2),
    color = dark_turquoise, key = "f = 1",  implicit (f = 1, x,0,2,y, 0,1),
    color = purple, key = "f = 1.5",  implicit (f = 1.5, x,0,2,y,0,1),
    color = red, key = "f = 1.8",  implicit (f = 1.8, x,0,2,y,0,2),
    color = green, key = "f = 2", points ([ [2,0] ])), wxplot_size = [680,680];    

We see from the plot that the maximum of f that satisfies x^2 + y^2 <= 4, x >=0, y>=0,
is the corner solution x = 2, y = 0.

gridSearch (f, [x,1.7, 2.3, 0.01], [y, -0.3, 0.3, 0.01], max, g);

KKT (f, max, g );

Check the first derivatives of the Lagrangian:

L : f + lam*g;
[Lx, Ly, Llam] : [diff(L,x), diff(L,y), diff(L,lam)];
subst ( [x = 2, y = 0, lam = 1/4], [f, x^2 + y^2, Lx, Ly, g] );

The optimal solution returned from KKT lies on the constraining border (the constraint is
binding) since g = 0.

## Dowling Prob. 13.39; Maximize Utility, 2 Constraints

Maximize the utility function u = x*y subject to the two constraints:
 1. budget constraint: 3 x + 4 y <= 144 --> 0 <= 144 - 3 x - 4 y, --> 144 - 3 x - 4 y >= 0,
 2. dietary constraint:  5 x + 2 y <= 120  -->  120 - 5 x - 2 y >= 0,
as well as x >= 0, y >= 0.

To use KKT(u, max, g1, g2) we need g1 >= 0, g2 >= 0.

Make a plot of the two constraint borders, and show u = constant curves (indifference
curves).

u : x*y;

wxdraw2d (yrange = [0,80], xrange = [0,80], line_width = 2, proportional_axes = 'xy,
    xlabel = "x", ylabel = "y", /* key = " feasible set ", 
    filled_func = true, fill_color = light_gray, */
    color = black, key = " need 3*x + 4*y <= 144", explicit (36 - 3*x/4, x, 0, 48),
    color = blue,  key = " need 5*x + 2*y <= 120",   explicit (60 - 5*x/2, x, 0, 24),
/*    filled_func = false, */
    color = dark_turquoise, key = "u = 100",  implicit (u = 100, x,0,30,y,0,30),
    color = red, key = "u = 200",  implicit (u = 200, x,0,30,y,0,30),
    color = purple, key = "u = 300",  implicit (u = 300, x,0,30,y,0,30))$
    

We see from the plot that the optimum (maximum) value of u = x*y is found near the point
where the constraint lines cross.

soln : solve ([3*x + 4*y = 144, 5*x + 2*y = 120], [x, y] )[1];

float(soln);

float (at (u, soln));

### gridSearch (f, xR, yR, maxmin, g1, g2 )

Use gridSearch (f(x,y), [x, x1, x2, dx], [y, y1, y2, dy], maxmin, g1, g2 )
in which if maxmin = max, then g1(x,y) >= 0 and g2(x,y) >= 0
and if maxmin = min, then g1(x,y) <= 0 and g2(x,y) <= 0.

g1 : 144 - 3*x - 4*y;
g2 : 120 - 5*x - 2*y;

gridSearch (u, [x, 13.5, 13.9, 0.01], [y, 25.4, 30, 0.01], max, g1, g2 );

We expect KKT to find this optimal solution with g1 >= 0, g2  >= 0.

KKT (u, max, g1, g2 );

Both constraints are binding (g1 = 0, g2 = 0) implies the optimum solution lies at the
intersection of these boundary curves.

L : u + lam1*g1 + lam2*g2;
[Lx, Ly, L1, L2] : [diff (L,x), diff(L,y),diff(L,lam1), diff (L, lam2)];
subst ( [x = 96/7, y = 180/7, lam1 = 60/49,lam2 = 216/49], [u, Lx, Ly, L1, L2, g1, g2] );

float(%);

## Chiang/Wainwright Ex. 1, Maximize Utility, 2 Constraints 

Chiang & Wainwright, p. 406-408

Maximize the utility u = x*y subject to the budget constraint x + y <= 100 and
a "ration constraint"   x <= 40.

Let g1 = 100 - x - y >= 0, and g2 = 40 - x >= 0.

In this problem with a fixed upper limit on good x, it is useful to make a plot showing the
region of feasible solution points (x,y), the kink in the border of the feasible set,
and a sample indifference curve u(x,y) = u0.

g1 : 100 - x - y;
g2 : 40 - x;
node1 : [x=0,y=100];
node2 : solve ([g1,g2], [x,y])[1];
node3 : [x = 40, y = 0];
y1 : rhs (solve(g1,y)[1]);

u(x,y) := x*y;

[u1,u2,u3] : map (lambda ([xx],at( u(x,y), xx)), [node1, node2, node3] );

wxdraw2d(line_width = 3, xrange = [0,120],yrange = [0,120],  key = "x + y = 100",
    /* two bounding curves */
    explicit (y1, x, 0, 100), color = black, key = "x = 40",
    parametric (40, t, t, 0, 120),
    filled_func = true, fill_color = light_gray, key = "feasible set",line_width = 1,
    explicit (y1,x,0,40), 
    /* three nodes */
    point_type = filled_circle, color = green,key = sconcat (node1," u = ",u1),  
    points ([ map('rhs, node1) ]),
     key = sconcat (node2," u = ",u2), color = magenta,points ([ map('rhs, node2) ]),
     key = sconcat (node3," u = ",u3), color = purple,points ([ map('rhs, node3) ]), 
    /* one u = constant curve */
    color = red,key = "u(x,y) = 2000",filled_func = false, line_width = 2, 
    explicit (2000/x, x, 18, 100))$

gridSearch (u,[x, 39,41,0.01],[y,59,61,0.01], max, g1, g2);

KKT (u, max, g1, g2 );

soln : [x=40,y=60,lam1=40,lam2=20];

L : u + lam1*g1 + lam2*g2;
[Lx, Ly, L1, L2] : [diff (L, x), diff (L, y), diff (L, lam1), diff (L, lam2) ];
subst (soln, [ u,x+y, x, Lx, Ly,L1, L2, g1, g2]);

The optimal solution has x = 40, x + y = 100, g1 = g2 = 0, u = 2400.

## Chiang/Wainwright Ex. 2, Minimize Cost, 2 Constraints    

p. 410 - 411.

Minimize C = (x - 4)^2 + (y - 4)^2 subject to
1.     2 x + 3 y >= 6 , 0 >= 6 - 2 x - 3 y
    g1(x,y) = 6 - 2 x - 3 y <= 0
2.   - 3 x - 2 y >= -12, 0 >= -12 + 3 x + 2 y,
    g2 (x, y) = -12 + 3 x + 2 y < = 0
Constraint #2 is equivalent to 3 x + 2 y <= 12, since if we multiply thru by (-1)
we muct reverse the direction of the inequality. (If -5 > -12 then 5 < 12)

The cost function is a circle of radius sqrt(C) with center at (x=4,y=4). With a protractor,
ruler and a piece of graph paper, draw the lines 2*x + 3*y = 6, 3*x + 2*y = 12, and 
circles of various radii centered at (4,4); the cost is the radius squared. You will quickly see
that the minimum cost corresponds to a point on the line 3*x + 2*y = 12, with a cost of
about 5.

C : (x - 4)^2 + (y - 4)^2;
g1 : 6 - 2*x - 3*y;
g2 : -12  + 3*x + 2*y;

gridSearch (C,[x,2,2.4,0.01],[y,2.6,3,0.01],min, g1, g2);

wxdraw2d (yrange = [0,8], xrange = [0,6], xlabel = "x", ylabel = "y",
   proportional_axes = 'xy, filled_func = true, fill_color  = light_gray,
 color = green, key = "need 3*x + 2*y <= 12", explicit (6 - 3*x/2,x,0,4),
 fill_color  = white, key = "need 2*x + 3*y >= 6", explicit (2 - 2*x/3,x,0,3),
color = dark_turquoise, key = "C = 7",  implicit (C = 7, x,1,4,y,1,5),
    color = red, key = "C = 6",  implicit (C = 6, x,1,4,y,1,5),
    color = purple, key = "C = 5",  implicit (C = 5, x,1,4,y,1,5)),wxplot_size = [680,680];

KKT (C, min, g1, g2 );

soln : [x=28/13, y=36/13, lam1=0, lam2=16/13];
L : C + lam1*g1 + lam2*g2;
[Lx : diff (L, x), Ly : diff (L,y), L1 : diff (L, lam1), L2 : diff (L, lam2) ];
float (subst (soln, [ C, Lx, Ly,L1, L2, g1, g2]));

The first constraint was g1 <= 0, or 2 x + 3 y >= 6. With g1 = -6.6, this is not binding:

at (2*x + 3*y,soln), numer;

The second constraint was g2 <= 0 or - 3 x - 2 y >= -12, and with g2 = 0, is a binding
constraint.

at (-3*x - 2*y, soln), numer;

## Maximize, 2 Constraints, Baldani, Bradfield, Turner

This example is from Baldani, Bradfield, and Turner, Mathematical Economics, 2nd. ed,
P. 291 - 293.

Maximize f(x,y) = x + y subject to x >= 0, y >= 0, x + 2 y <= 5, 2 x + y <= 6.

Define g1(x,y) and g2(x,y) so g1 >= 0, g2 >= 0.

f : x + y;
g1 : 5 - x - 2*y;
g2 : 6 - 2*x - y;

Where do the lines g1 = 0, g2 = 0 cross?

cross : solve ([g1,g2],[x,y])[1];

float(cross);

float (at(f, cross));

Graphical exploration; show feasible (x,y) set in light gray, and some
curves f = constant.

wxdraw2d(line_width = 2, xrange = [0,5], yrange = [0,6],   xlabel = "x", ylabel = "y",
    proportional_axes = 'xy,   key = "x + 2*y = 5", explicit (5/2 - x/2,x, 0, 5), 
    color = green, key = "2*x + y = 6",
    explicit (6 - 2*x, x, 0, 3), key = "feasible set",
    filled_func = true, fill_color = gray,
    explicit (if x < 7/3 then 5/2 - x/2 else 6 - 2*x, x, 0, 3),
    color = black, key = "f = 2",  implicit (f = 2, x,0,4,y,0,4),
    color = red, key = "f = 3",  implicit (f = 3, x,0,4,y,0,4),
    color = purple, key = "f = 4",  implicit (f = 4, x,0,4,y,0,4)),wxplot_size = [680,680];
    

gridSearch (f,[x,2.1,2.4,0.01],[y,1.1,1.5,0.01], max,  g1, g2 );

KKT (f, max, g1, g2 );