% PC profit example
%
% This code solves the PC profit optimization problem we covered in class,
% and calculates a few sensitivity coefficients.
%
% You can either run the code all at once, using the 'run' command, or in
% sections, using 'run section'. Sections are created using '%%'.
                                                
% Note: if you ever want to start the file over, type 'clear' to delete the
% workspace. 'clc' clears the command window.

%% Define variables and parameters/constants

% Note: the semicolon at the end of each line suppresses output.

n = 10000;  % Current sales [units/month]
m = 700;    % Manufactoring cost [$/unit]
w = 950;    % Wholesale unit price [$/unit]
d = 100;    % Amount of test discount [$/unit]
r = 0.5;    % Percent sales increase based on test discount [-]
a = 50000;  % Current advertizing budget [$/month]
b = 10000;  % Amount of test ad budget increase [$/month]
c = 200;    % Increase in sales due to ad budget test [units/month]
A = 100000; % Cap on ad budget [$/month]

%% Plot the model

D = linspace(0,3,200);      % Create a test vector for D and B, since we haven't specified it. 
B = linspace(-100,100,200);
% D = linspace(-10,10,200);     % This range produces a better-looking plot of the surface P
% B = linspace(-1000,1000,200);
[Dmesh, Bmesh] = meshgrid(D,B); % Combine D and B into a grid that we can plot
N = n*(1 + r*Dmesh) + c*Bmesh;
q = w - m - d*Dmesh;
P = N.*q - a - b*Bmesh;

mesh(D,B,P)
xlabel('D'), ylabel('B'), zlabel('P')

%% Symbolic differentiation

% Create symbolic variables and equations
syms sr sD sB       % This initializes sr and sD as symbolic variables
sN = n*(1 + sr.*sD) + c*sB;
sq = w - m - d*sD;
sP = sN.*sq - a - b*sB; 

% Differentiate: 
dPdD = diff(sP,sD)  % d(total profit)/d(discount)
simplify(dPdD)
dPdB = diff(sP,sB)  % d(total profit)/d(budget)
simplify(dPdB)

%% Find optimal discount

% Set dPdD=0 and dPdB=0 and solve for sD and sB
% This gives us the optimal discount and ad budget as a function of sr.
[sD, sB] = solve([dPdD==0,dPdB==0],[sD,sB])

%% Lagrange multipliers

%% Maximize on the boundary B=10
clear sD sB
syms sD sB lam1

% Solve the system dPdD=0, dPdB=lambda, B=10 for lambda
[sD, sB, lam1] = solve([dPdD==0,dPdB==lam1,sB==10],[sD,sB,lam1])

sD1 = subs(sD,sr,0.5)
lam1 = subs(lam1,sr,0.5)

%% Maximize on the boundary B=0
clear sD sB
syms sD sB lam2

% Solve the system dPdD=0, dPdB=lambda, B=10 for lambda
[sD, sB, lam2] = solve([dPdD==0,dPdB==lam2,sB==0],[sD,sB,lam2])

sD2 = subs(sD,sr,0.5)
lam2 = subs(lam2,sr,0.5)

%% Plot level sets of P(D,B)

D = linspace(-1,1,200);      % Create a test vector for D and B, since we haven't specified it. 
B = linspace(-1,11,200);
[Dmesh, Bmesh] = meshgrid(D,B); % Combine D and B into a grid that we can plot
N = n*(1 + r*Dmesh) + c*Bmesh;
q = w - m - d*Dmesh;
P = N.*q - a - b*Bmesh;

figure, hold on
set(gca,'FontSize',18)
[C,h] = contourf(D,B,P);
contour(D,B,Bmesh,[0,0],'Color','k')
contour(D,B,Bmesh,[10,10],'Color','k')
xlabel('D'), ylabel('B')
clabel(C,h)