%% Linear Static Models
%
% Linear programming (LP) is a tool for solving optimization problems. In 1947, George Dantzig developed an efﬁcient method, the simplex algorithm, for solving linear programming problems (also called LP). Since the development of the simplex algorithm, LP has been used to solve optimization problems in industries as diverse as banking, education, forestry, petroleum, and trucking. In a survey of Fortune 500 ﬁrms, 85% of the respondents said they had used linear programming.
%
%% Example 1: Giapetto’s Woodcarving Problem
% This program solves Giapetto's Woodcarving problem with the matlab function "linprog".
%% Problem Statement:
% Giapetto’s Woodcarving, Inc., manufactures two types of wooden toys: soldiers and trains. A soldier sells for $27 and uses $10 worth of raw materials. Each soldier that is manufactured increases Giapetto’s variable labor and overhead costs by $14. A train sells for $21 and uses $9 worth of raw materials. Each train built increases Giapetto’s variable labor and overhead costs by $10. The manufacture of wooden soldiers and trains requires two types of skilled labor: carpentry and ﬁnishing. A soldier requires 2 hours of ﬁnishing labor and 1 hour of carpentry labor. A train requires 1 hour of ﬁnishing and 1 hour of carpentry labor. Each week, Giapetto can obtain all the needed raw material but only 100 ﬁnishing hours and 80 carpentry hours. Demand for trains is unlimited, but at most 40 soldiers are bought each week. Giapetto wants to maximize weekly proﬁt (profit = revenues - costs). Formulate a mathematical model of Giapetto’s situation that can be used to maximize Giapetto’s weekly proﬁt.
%% Solution:
%  In developing the Giapetto model, we explore characteristics shared by all linear programming problems.
%% Decision Variables
%  We begin by deﬁning the relevant decision variables. In any linear programming model, the decision variables should completely describe the decisions to be made (in this case, by Giapetto). Clearly, Giapetto must decide how many soldiers and trains should be manufactured each week. With this in mind, we deﬁne,
% x_1 = number of soldiers produced each week
% x_2 = number of trains produced each week
%% Objective Function
%  In any linear programming problem, the decision maker wants to maximize (usually revenue or profit) or minimize (usually costs) some function of the decision variables. The function to be maximized or minimized is called the objective function. For the Giapetto problem, we note that fixed costs (such as rent and insurance) do not depend on the values of x_1 and x_2. Thus, Giapetto can concentrate on maximizing (weekly revenues) - (raw material purchase costs) - (other variable costs).
% Giapetto’s weekly revenues and costs can be expressed in terms of the decision variables x_1 and x_2. It would be foolish for Giapetto to manufacture more soldiers than can be sold, so we assume that all toys produced will be sold. Then
% Weekly revenues = weekly revenues from soldiers +
%                                 weekly revenues from trains
%                                = \Big(\frac{dollars}{soliders}\Big)\Big(\frac{sloiders}{week}\Big)+ \Big(\frac{dollars}{train}\Big)\Big(\frac{trains}{week}\Big)
%                                = 27x_1+21x_2
% Also,
% Weekly raw material costs   = 10x_1+9x_2
% Other weelky variable costs = 14x_1+10x_2
% Then Giapetto wants to maximize,
% max \ \pi = (27x_1+21x_2)-(10x_1+9x_2)-(14x_1+10x_2)=3x_1+2x_2
% Another way to see that Giapetto wants to maximize 3x_1+2x_2 is to note that,
% Weekly revenues = weekly contribution to profit from soldiers - weekly non-fixed costs +
%                                 weekly contribution to profit from trains
%                             = \Big(\frac{contribution \ to \ profit}{soliders}\Big)\Big( \frac{soliders}{week} \Big) + \Big(\frac{contribution \ to \ profit}{trains}\Big)\Big(\frac{trains}{week}\Big)
% Also,
% Contribution to Profit = 27 - 10 - 14 = 3
% Contribution to profit = 21 - 9 - 10 = 2
% Then as before we obtain,
% Weekly revenues - weekly non-fixed costs = 3x_1+2x_2
% Thus, Giapetto’s objective is to choose x_1 and x_2 to maximize 3x_1+2x_2. We use the variable z to denote the objective function value of any LP. Giapetto’s objective function is
% Maximize \ z =3x_1 + 2x_2 =\pi \quad \quad\quad \quad (1)
% (In the future, we will abbreviate “maximize” by max and “minimize” by min.) The coefficient of a variable in the objective function is called the objective function coefficient of the variable. For example, the objective function coefficient for x_1 is 3, and the objective function coefficient for x_2 is 2. In this example (and in many other problems), the objective function coefficient for each variable is simply the contribution of the variable to the company’s profit.
%% Constraints
%  As x_1 and x_2 increase, Giapetto’s objective function grows larger. This means that if Giapetto were free to choose any values for x_1 and x_2, the company could make an arbitrarily large profit by choosing x_1 and x_2to be very large. Unfortunately, the values of x_1 and x_2 are limited by the following three restrictions (often called constraints):
% Constraint 1: Each week, no more than 100 hours of finishing time may be used.
% Constraint 2: Each week, no more than 80 hours of carpentry time may be used.
% Constraint 3: Because of limited demand, at most 40 soldiers should be produced each week.
% The amount of raw material available is assumed to be unlimited, so no restrictions have been placed on this.
% The next step in formulating a mathematical model of the Giapetto problem is to express Constraints 1–3 in terms of the decision variables x_1 and x_2. To express Constraint 1 in terms of x_1 and x_2, note that,
%  Total \ finishing \ hours = \Big(\frac{finishing \ hours}{soliders}\Big)\Big(\frac{soliders \ made}{week}\Big) + \Big(\frac{finishing \ hours}{trains}\Big)\Big(\frac{trains \ made}{week}\Big)
%                                      = 2x_1+x_2
% Now Constrain 1 may be expressed by
% 2x_1+x_2\leq100 \quad\quad\quad\quad\quad\quad(2)
% Note that the units of each term in (2) are finishing hours per week. For a constraint to be reasonable, all terms in the constraint must have the same units. Otherwise one is adding apples and oranges, and the constraint won’t have any meaning.
% To express Constraint 2 in terms of x_1 and x_2, note that
% Total \ carpentry \ hours = \Big(\frac{carpentry \ hours}{soliders}\Big)\Big(\frac{soliders}{week}\Big) + \Big(\frac{carpentry}{trains}\Big)\Big(\frac{trains}{week}\Big)
%                                       = x_1+x_2
% Then Constraint 2 may be written as
% x_1+x_2\leq80 \quad\quad\quad\quad\quad\quad(3)
% Again, note that the units of each term in (3) are the same (in this case, carpentry hours per week).
% Finally, we express the fact that at most 40 soldiers per week can be sold by limiting the weekly production of soldiers to at most 40 soldiers. This yields the following constraint:
% x_1\leq40 \quad\quad\quad\quad\quad\quad(4)
% Thus (2)–(4) express Constraints 1–3 in terms of the decision variables; they are called the constraints for the Giapetto linear programming problem. The coefficients of the decision variables in the constraints are called technological coefficients. This is because the technological coefficients often reflect the technology used to produce different products. For example, the technological coefficient of x_2 in (3) is 1, indicating that a soldier requires 1 carpentry hour. The number on the right-hand side of each constraint is called the constraint’s right-hand side (or rhs). Often the rhs of a constraint represents the quantity of a resource that is available.
%% Sign Restrictions
% To complete the formulation of a linear programming problem, the following question must be answered for each decision variable: Can the decision variable only assume nonnegative values, or is the decision variable allowed to assume both positive and negative values?
% If a decision variable x_i can only assume nonnegative values, then we add the sign restriction x_i\geq0. If a variable x_i can assume both positive and negative (or zero) values, then we say that x_i is unrestricted in sign (often abbreviated urs). For the Giapetto problem, it is clear that x_1\geq0 and x_2\geq0. In other problems, however, some variables may be urs. For example, if x_i represented a firm’s cash balance, then x_i could be considered negative if the firm owed more money than it had on hand. In this case, it would be appropriate to classify x_i as urs.
% Combining the sign restrictions x_1\geq0 and x_2\geq0 with the objective function (1) and Constraints (2)–(4) yields the following optimization model:
% max \ z = 3x_1+2x_2 \quad\quad\quad\quad(Objective \ function)
% subject to (s.t)
%                                                 2x_1+x_2\leq100 \quad\quad\quad\quad(Finishing \ constraints)
%                                                 x_1+x_2\leq80 \quad\quad\quad\quad\quad\quad(Finishing \ constraint)
%                                                 x_1\quad\quad\leq40 \quad\quad\quad\quad\quad\quad(Constraint \ on \ demand \ for \ soliders)
%                                                 x_1\quad\quad\geq0 \quad\quad\quad\quad\quad\quad(sign \ restriction)
%                                                         x_2\geq0 \quad\quad\quad\quad\quad\quad(sign \ restriction)
% “Subject to” (s.t.) means that the values of the decision variables x_1 and x_2 must satisfy all constraints and all sign restrictions.
%
%% DEFINITION
% A linear programming problem (LP) is an optimization problem for which we do the following:
% We attempt to maximize (or minimize) a linear function of the decision variables. The function that is to be maximized or minimized is called the objective function.
% The values of the decision variables must satisfy a set of constraints. Each constraint must be a linear equation or linear inequality.
% A sign restriction is associated with each variable. For any variable x_i, the sign restriction specifies that x_i must be either nonnegative (x_i\geq0) or unrestricted in sign (urs).
% Because Giapetto’s objective function is a linear function of x_1 and x_2, and all of Giapetto’s constraints are linear inequalities, the Giapetto problem is a linear programming problem. Note that the Giapetto problem is typical of a wide class of linear programming problems in which a decision maker’s goal is to maximize profit subject to limited resources.
%
%

% This program solves Giapetto's Woodcarving problem with the matlab
% function "linprog."

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the parameter vector in the objective function.

clear;clc;

f = [-3 -2]; % profits per unit of output

% linprog by defaut is a minimization algoritm (type: 'help linprog' into the command 
% window to see more details). We can change this to a maximization algorithm by changing 
% the sign of the function 'f' or setting the signs of the paramaters of the function 
% to thier opposite.
 
% Define matrix A containing technological paremeters in the ineqaulity
% constraints.

A = [2 1;
     1 1;
     1 0;
    -1 0;
    0 -1];
%  
% Define vector b.

b = [100; 80; 40; 0; 0];

% Call a linear programming routine.

[x, fval, exitflag, output, lambda] = linprog(f, A, b);
x

%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The problem can be set up in a different way as follows

% Define the parameter vector in the objective function.

clc;clear;

f = [-3 -2];
 
% Define matrix A containing technological paremeters in the ineqaulity
% constraints.

A = [2 1;
     1 1];
%  
% Define vector b.

b = [100; 80];

% Define the lower and upper bound vectors

lb = [0; 0];
ub = [40; Inf];


% Call a linear programming routine.

[x, fval, exitflag, output, lambda] = linprog(f,A,b,[],[],lb,ub);
x

%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we can change the constraints

% Define the parameter vector in the objective function.
clc;clear;

f = [-3 -2];
 
% Define matrix A containing technological paremeters in the ineqaulity
% constraints.

A = [2 1;
     1 1];
%  
% Define vector b.

b = [101; 80];

% Define the lower and upper bound vectors

lb = [0; 0];
ub = [40; Inf];


% Call a linear programming routine.

[x, fval, exitflag, output, lambda] = linprog(f,A,b,[],[],lb,ub);
x

%
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program determines Giapetto's production choices and weekly profit
% when the price of wooden soldiers varies in the rage of $20 - 40$.

% To run the following script, you need to save the m-file 'giap' in the
% same folder in which the present script is saved.

clc;clear;

P = 20:40;              % Define the range of the price of wooden soldiers.
z = zeros (8,numel(P)); % Create a zero matrix for storing results.

for k = 1:numel(P)
P1 = P(k);   % price of wooden soldiers
Cm1 = 10;    % Material costs per wooden soldier
Cvf1 = 14;   % Variable costs and overhead per wooden soldier
Lf1 = 2;     % Finishing labor requirements per wooden soldier
Lc1 = 1;     % Carpentry labor requirements per wooden soldier

P2 = 21;     % Price of wooden trains
Cm2 = 9;     % Material cost per wooden train
Cvf2 = 10;   % Variable costs and overhead per wooden train
Lf2 = 1;     % Finishing labor requirements per wooden train
Lc2 = 1;     % Carpentry labor requirements per wooden train

Lfmax = 100; % Maximium finishing labor
Lcmax = 80;  % Maximum carpentry labor
D1max = 40;  % Maximum weekly demand for wooden soldiers

% Call the function, giap, to solve for the optimal outputs of woodern 
% soldiers and trains.  Save optimal outputs, weekly profit, and Lagrangian
% multiplies in output matrix z.

[x, fval, exitflag, lambda] = giap(P1, Cm1, Cvf1, Lf1, Lc1, P2, Cm2, Cvf2, Lf2, Lc2,...
    Lfmax, Lcmax, D1max);

z ([1 2], k) = x;
z (3, k) = -fval;
z (4, k) = lambda.ineqlin(1);
z (5, k) = lambda.ineqlin(2);
z (6, k) = lambda.upper(1);
z (7, k) = P(k);
z (8, k) = exitflag;

end

rowNames = {'x1','x2','fval','lambda_ineqlin_1','lambda_inequal_2','lambda_upper',...
    'price_soldiers', 'exitflag'};
results = array2table(z,'RowNames',rowNames);

% Plotting the solution functions and value function with respect to price
% of wooden soldiers.

%subplot(3,2,1)
plot(P, z(1,:));
ylabel('Weekly Output of Wooden Soldiers')
xlabel('Price of Woodern Soldiers')

%subplot(3,2,2)
plot(P, z(2,:));
ylabel('Weekly Output of Wooden Trains')
xlabel('Price of Woodern Trains')

%subplot(3,2,3)
plot(P, z(3,:));
ylabel('Weekly Profit')
xlabel('Price of Woodern Soldiers')

%subplot(3,2,4)
plot(P, z(4,:));
ylabel('Shadow Price -- Finishing Labor')
xlabel('Price of Woodern Soldiers')

%subplot(3,2,5)
plot(P, z(5,:));
ylabel('Shadow Price -- Capertry Labor')
xlabel('Price of Woodern Soldiers')

%subplot(3,2,6)
plot(P, z(6,:));
ylabel('Shadow Price -- Demand Constraint - Woodern Soldiers')
xlabel('Price of Woodern Soldiers')

%
%% Example 1 - Extension: The Dual Approach - Giapetto’s Woodcarving Problem
%
% In this first part of the Giapetto Woodcarving problem we formulated a situation where we desired to maximize profits subject to a set of constraints. However, as you recall from your microeconomics we can also formulate the dual problem where we minimize costs subject to a set of constrains.
% Primal Problem:
% \[\begin{array}{r}
% \max \;z = \sum\limits_{j = 1}^n {{c_j}{x_j}\quad s.t.\quad \sum\limits_{j = 1}^n {{a_{ij}}{x_j} \le {b_j}\quad \forall \;i = 1,2,...,m} } \\
% {x_j} \ge 0\quad \forall \;j = 1,2,...,n
% \end{array}\]
%
% Dual Problem:
% \[\begin{array}{r}
% \min \;v = \sum\limits_{i = 1}^m {{b_i}{y_i}\quad s.t.\quad \sum\limits_{i = 1}^m {{a_{ij}}{y_i} \ge {c_j}\quad \forall \;i = 1,2,...,m} } \\
% {y_i} \ge 0\quad \forall \;j = 1,2,...,n
% \end{array}\]
%
% Using the same information given in the above example, see if you can reformulate the primal problem (profit max.) to the dual problem (cost min.):
%% Problem Statement:
% Giapetto’s Woodcarving, Inc., manufactures two types of wooden toys: soldiers and trains. A soldier sells for $27 and uses $10 worth of raw materials. Each soldier that is manufactured increases Giapetto’s variable labor and overhead costs by $14. A train sells for $21 and uses $9 worth of raw materials. Each train built increases Giapetto’s variable labor and overhead costs by $10. The manufacture of wooden soldiers and trains requires two types of skilled labor: carpentry and ﬁnishing. A soldier requires 2 hours of ﬁnishing labor and 1 hour of carpentry labor. A train requires 1 hour of ﬁnishing and 1 hour of carpentry labor. Each week, Giapetto can obtain all the needed raw material but only 100 ﬁnishing hours and 80 carpentry hours. Demand for trains is unlimited, but at most 40 soldiers are bought each week.
% Hint: Think of the objective function as the "cost of capacity" so our objective function will be given as,
% min \ v = 100y_1+80y_2
% If the problem is formulated correctly you should see a relationship between the optimal values of y_i and the lambda values of the primal problem. In addition, you should see a relationship between the optimal values of x_i from the primal problem and the lambda values from the dual approach. This is not an accidental relationship.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The Dual problem can be set up as follows
clc;clear;

% Define the parameter vector in the objective function.

f = [100 80];
 
% Define matrix A containing technological paremeters in the ineqaulity
% constraints.

A = [];
%  
% Define vector b.

b = [];

% Define the lower and upper bound vectors

lb = [];
ub = [];


% Call a linear programming routine.

[x, fval, exitflag, output, lambda] = linprog(f,A,b,[],[],lb,ub);
x

%
%% Shadow Price (lambda or lagrangian multiplier)
% A shadow price represents the maximal price at which we should be willing to buy additional units of a resource. It also represents the minimal price at which we should be willing to sell units of the resource. A shadow price might therefore be thought of as the value per unit of a resource (marginal). Remarkably, if we compute the value of our entire stock of resources based on shadow prices, we get our optimal profit! For instance, in our example of building soldiers and trains, we have
% 1(100) + 1(80) = 180
% It is not a coincidence that the minimal cost in the dual problem equals the optimal profit in the primal and that the optimal solution of the dual is the vector of shadow prices – these are fundamental relations between the primal and the dual.
%%
%% References:
%
% Winston, W. L., & Goldberg, J. B. (2004). Operations research: applications and algorithms. Duxbury Press.
%