Constrained Whales (Prb 2.4.2)

Variables: General Case

x=number of blue whales
y=number of fin whales

Gx=birth rate of blue whales in whales per year
Gy=birth rate of fin whales in whales per year
Parameters:

r1=natural growth rate of blue whales (blue whales per year per blue whale)
r2= natural growth rate of fin whales(fin whales per year per fin whale)

K1=environmental carrying capacity for blue whales (in blue whales)
K2=environmental carrying capacity for fin whales(in fin whales)

alpha=competition coefficient, assumed to be the same for blues and fins

Values of parameters at which we will compute various quantities:

r1=.05
r2=.08
K1=150000
K2=400000
alpha=10^(-8)

Relations: These are general; later we substitute in the parameter values as appropriate.

> restart;
Gx:=r1*x*(1-x/K1)-alpha*x*y;
Gy:=r2*y*(1-y/K2)-alpha*x*y;

Objective: Determine the population levels x,y which maximize subject to the conditions .

The objective function here is x+y, the decision variables are x and y, and the gradient never vanishes. Hence we will find the maximum point on the boundary of the region. We have to find the region. Here are some macros for later use.

> Cons:={r1=5/100,r2=8/100, K1=150000, K2=400000,alpha=10^(-8)}; #to be used to fix the parameters
numcons:=a->subs(Cons,a); #fixes the parameters as specified by Cons;
sens:=(f,x)->diff(f,x)*x/f; #sens(f,x) is the sensitivity of f wrt x
lm:=(f,g)->solve({diff(f,x)=lambda*diff(g,x),diff(f,y)=lambda*diff(g,y),g},{x,y,lambda});
# in case we need to find some lagrange multipliers.

We plot the feasible region. We know what looks like, and have to find what and look like. These are simple questions. E.g., the level set Gx=0 separates regions where 0<Gx and Gx>0, so we figure this out. Gx=0 requires either x=0 or (1-x/K1)-alpha*y=0, and this is a line. We plot these lines arising from Gx=0 and Gy=0 together with a level set of x+y.

> plots[implicitplot]({numcons(Gx/x),numcons(Gy/y),x+y=300000},x=0 .. numcons(K1),y=0..numcons(K2),labels=['x','y']);

The top line is part of Gy=0 since it contains the point x=0, y=K2. We see that the linear objective function is not parallel to the bounding lines, so its gradient is nowhere parallel to the gradients of any bounding line x=0, y=0, Gx=0, Gy=0 on the boundary. We need only check the corners to find the maximum. Let us confirm via a Lagrange multiplier check:

> lm(x+y,numcons(Gx/x));
lm(x+y,x);lm(x+y,y);
lm(x+y,numcons(Gy/y));

There are no solutions. We compute the values of x+y at the corners:

> numcons(subs({x=0,y=K2},x+y));
numcons(subs({x=K1,y=0},x+y));
upper:=solve({Gx/x,Gy/y},{x,y}); # the coordinates of the upper corner
evalf(numcons(subs(upper,x+y)));

Thus the maximum point and the maximum value are
evalf(numcons(upper));
evalf(numcons(subs(upper,x+y)));

It is clear that the situation will be basically the same for nearby values of the environment variables, . The point "upper" gives the maximum point and

> maxval:=subs(upper,x+y);

gives the maximum value. It is routine to compute the sensitivities of maxval to the parameters and we do it as a display of power:

> assign(upper);

> evalf(numcons(sens(x,r1)));
evalf(numcons(sens(x,r2)));
evalf(numcons(sens(y,r1)));
evalf(numcons(sens(y,r2)));
evalf(numcons(sens(x+y,r1)));
evalf(numcons(sens(x+y,r2)));
evalf(numcons(sens(x,K1)));
evalf(numcons(sens(x,K2)));
evalf(numcons(sens(y,K1)));
evalf(numcons(sens(y,K2)));
evalf(numcons(sens(x+y,K1)));
evalf(numcons(sens(x+y,K2)));

> unassign('x','y');

To see what happens as a function of alpha, we animate the situation. First, let us solve the various equations we plot for y, leaving alpha variable:, incorporating an animation parameter t:

> Cons:={r1=5/100,r2=8/100, K1=150000, K2=400000,alpha=t*10^(-8)}:
y1:=solve(numcons(Gx/x),y);
y2:=solve(numcons(Gy/y),y);
y3:=300000-x;

> with(plots):

> animate({y1,y2,y3},x=0 .. numcons(K1),t=1..13,
color=cyan,view=[0..numcons(K1),0...numcons(K2)]);

We conclude that as alpha increases there is a first alpha such that the line Gx/x=0 intersects the y axis at the same point as Gy/y=0 (y=400000) and at that point, the optimal solution is x=0, y=400000. We calculate the alpha:

> Cons:={r1=5/100,r2=8/100, K1=150000, K2=400000,alpha=alpha};

> solve(subs({x=0,y=400000},simplify(numcons(Gx/x))),alpha);

> evalf(%);

>

The phenomenon occurs at alpha=.125*10^(-6).