function prodDynamic

%The horizon to optimize over:
horizon = ????; 

numStates = 12;

initStatesFixed = [1 2 3 4 5 6 7];
x0 = zeros(length(initStatesFixed), 1);

%Make the cost function:
p = [21 18 -0.5 -0.5];
c = costVector(p, horizon, numStates);

%Make the constraints:
[A b Aeq beq LB UB] = constraints(x0, initStatesFixed, horizon, numStates);

%Solve the optimization problem:
x = linprog(-c, A, b, Aeq, beq, LB, UB);
%Arrange the solution so that the t-th column in X contains the
%states for time t:
X = reshape(x, numStates, horizon)


%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = costVector(p, horizon, numStates)
c = zeros(horizon*numStates, 1);
c(1+(0:horizon-1)*numStates) = p(1)*ones(horizon,1); 
c(2+(0:horizon-1)*numStates) = p(2)*ones(horizon,1);
c(5+(0:horizon-1)*numStates) = -0.5*ones(horizon,1);
c(6+(0:horizon-1)*numStates) = -0.5*ones(horizon,1);
c(8+(0:horizon-1)*numStates) = -ones(horizon,1);
c(11+(0:horizon-1)*numStates) = p(3)*ones(horizon,1);
c(12+(0:horizon-1)*numStates) = p(4)*ones(horizon,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A b Aeq beq LB UB] = constraints(x0, iSF, horizon, numStates)
LB = zeros(horizon*numStates, 1);
LB(8+(0:horizon-1)*numStates) = 100*ones(horizon,1); %x_8(t) >= 100
LB(9+(0:horizon-1)*numStates) = -100*ones(horizon,1); %x_9(t) >= -100
LB(10+(0:horizon-1)*numStates) = -100*ones(horizon,1); %x_10(t) >= -100

UB = zeros(horizon*numStates, 1);
UB(1+(0:horizon-1)*numStates) = Inf*ones(horizon,1); %x_1(t) unbounded
UB(2+(0:horizon-1)*numStates) = Inf*ones(horizon,1); %x_2(t) unbounded
UB(3+(0:horizon-1)*numStates) = 500*ones(horizon,1); %x_3(t) <= 500
UB(4+(0:horizon-1)*numStates) = 250*ones(horizon,1); %x_4(t) <= 250
UB(5+(0:horizon-1)*numStates) = Inf*ones(horizon,1); %x_5(t) unbounded
UB(6+(0:horizon-1)*numStates) = 1000*ones(horizon,1); %x_6(t) <= 1000
UB(7+(0:horizon-1)*numStates) = 750*ones(horizon,1); %x_7(t) <= 750
UB(8+(0:horizon-1)*numStates) = Inf*ones(horizon,1); %x_8(t) unbounded
UB(9+(0:horizon-1)*numStates) = 100*ones(horizon,1); %x_9(t) <= 100
UB(10+(0:horizon-1)*numStates) = 100*ones(horizon,1); %x_10(t) <= 100
UB(11+(0:horizon-1)*numStates) = Inf*ones(horizon,1); %x_11(t) unbounded
UB(12+(0:horizon-1)*numStates) = Inf*ones(horizon,1); %x_12(t) unbounded

A = [];
b = [];
for t = 0:horizon-1,
  %7x_1 + 10x_2 - x_3 - x_4 - x_5 <= 0
  addA = zeros(1,horizon*numStates);
  addA(numStates*t+[1 2 3 4 5]) = [7 10 -1 -1 -1];
  A = [A; addA];
  b = [b; 0];
  
  %16x_1 + 12x_2 - x_6 - x_7 <= 0
  addA = zeros(1,horizon*numStates);
  addA(numStates*t+[1 2 6 7]) = [16 12 -1 -1];
  A = [A; addA];
  b = [b; 0];
  
  %1.5x_7 - 425 <= x_8
  addA = zeros(1,horizon*numStates);
  addA(numStates*t+[7 8]) = [1.5 -1];
  A = [A; addA];
  b = [b; 425];
  
  %x_9 <= x_11
  addA = zeros(1,horizon*numStates);
  addA(numStates*t+[9 11]) = [1 -1];
  A = [A; addA];
  b = [b; 0];
  
  %-x_9 <= x_11
  addA = zeros(1,horizon*numStates);
  addA(numStates*t+[9 11]) = [-1 1];
  A = [A; addA];
  b = [b; 0];
  
  %x_10 <= x_12
  addA = zeros(1,horizon*numStates);
  addA(numStates*t+[10 12]) = [1 -1];
  A = [A; addA];
  b = [b; 0];
  
  %-x_10 <= x_12
  addA = zeros(1,horizon*numStates);
  addA(numStates*t+[10 12]) = [-1 1];
  A = [A; addA];
  b = [b; 0];
end

%x(0) = x0:
Aeq = zeros(length(iSF), horizon*numStates);
Aeq(iSF, iSF) = eye(length(iSF));
beq = x0;

for t = 0:horizon-2,
  %x_3(t+1) = x_3(t) + x_9(t)
  addAeq = zeros(1,horizon*numStates);
  addAeq(numStates*t+[3+numStates 3 9]) = [1 -1 -1];
  Aeq = [Aeq; addAeq];
  beq = [beq; 0];
  
  %x_4(t+1) = x_4(t) + x_10(t)
  addAeq = zeros(1,horizon*numStates);
  addAeq(numStates*t+[4+numStates 4 10]) = [1 -1 -1];
  Aeq = [Aeq; addAeq];
  beq = [beq; 0];
end
