This is my original code which is only slightly different from the picture provided in initial post for comparison with mathematical model. This should work without any errors (except for being infeasible)
{
clear
close all
clc
%% Data/Parameters used in the Model
%%
%Number of time periods
P = 32;
%Number of states machines can be in, off, on, cooling, idle, processing
S = 5;
%Total number of machines in the job shop
M = 2;
%Total number of jobs to be processed in the time span
J = 4;
%Demand charge per kW
DemandCharge = [10];
%Number of states in which the transition must take a pre-defined period of time before it can transition into a different state
%Rolling window of time for peak load
RollingWindow = [1];
%Process time for each job, in time periods
ProcTime = [10 9 10 9];
%Transitions allowed from Off (1) state, can stay off or turn On:
%Transitions allowed from Idle(3) state, can stay idle, start processing, or cool down
%Transitions allowed from Processing(4) state, stay in processing, or idle, or cool down
%Transitions allowed from Cool Down(5) state, can stay in coold down or turn off
TA{1} = [1 2];
TA{2} = [2 3 4];
TA{3} = [3 4 5];
TA{4} = [3 4 5];
TA{5} = [1 5];
%State in which machine is allowed to process jobs, that is, processing state (SET)
ProcState = [4];
%States in which transition takes predefined time period to complete before it can transition to a different state (also SET)
TransitionStates = [2 5];
%States in which machine is turned off or will be turned off (SET)
OffState = [1 5];
%States in which machine is turned off or will be turned off (SET)
%Time period that must be elapsed in the particular state before the machine is allowed to transition to a different state
%If machine is in state 1(ON) it must stay ON before it can idle or start processing job
%Machine needs to elapse a time period of 1 before it will be turned off
TransitionTime = [0 5 0 0 1];
%Electricity price in respective time period
EnergyPrice = [0.034 0.027 0.032 0.035 0.035 0.037 0.036 0.039 0.039 0.040 0.052 0.074 0.050 0.094 0.100 0.084 0.090 0.074 0.094 0.084 0.027 0.091 0.045 0.048 0.045 0.040 0.031 0.039 0.070 0.027 0.080 0.092];
%Energy consumption kWh in respective state
EnergyDemand = [0 5 2 6 1];
%due date for jobs
DueDates = [20 12 15 32];
%early penalty
EarlyPenalty = [2 3 2 1];
%tardiness penalty for jobs
TardinessPenalty = [3 3 1 4];
%% Call cvx
cvx_begin
cvx_solver Mosek
%% Variables used in the model
%%
variable W(M,S,P) binary; %Machine i in power state s during period p
variable X(M,J,P) binary; %Machine i is processing job j during period p
variable Y(M,J,P) binary ; %Machine i starts to process job j during period p
variable Z(M,S,P) binary; %Machine i begins transition state a during period p, a must be a power state which allows transition from one state to other, such as transition from off to on is allowed but off to processing is not
variable d Nonnegative; %max peak energy in an Rw period interval
variable EarlinessJob(J,1) Nonnegative; %Earliness of job j
variable T(J,1) Nonnegative; %Lateness of job j
variable CompletionTime(J,1) Nonnegative;
variable Tardiness Nonnegative;
variable Earliness Nonnegative;
variable DemandCost Nonnegative;
variable MainCost Nonnegative;
%% Constraints and Objective
%% Objective function
minimize ((DemandCost) + (MainCost) + (Tardiness + Earliness));
subject to
%% Constraints
%%
%Objective Constraints
%Tardiness
for j = 1:J
Tardiness == sum(T(:) .* TardinessPenalty(:));
end
for j = 1:J
Earliness == sum(EarlinessJob(:) .* EarlyPenalty(:));
end
DemandCost == DemandCharge * d;
for i = 1:M
for s = 1:S
for p = 1:P
MainCost == sum(sum(sum(W(i,s,p) .* EnergyDemand(s) .* EnergyPrice(p))));
end
end
end
%constraint 1: Ensures Machine is in processing state while processing the job, the element k will be in Pr set to ensure that
for i = 1:M
for s = ProcState
for p = 1:P
sum(X(:,j,:)) == sum(W(:,ProcState,:));
end
end
end
%constraint 2: Machine is permitted to be in just one state during any period
for i = 1:M
for p = 1:P
for s = 1:S
sum (W(:,s,:)) == 1;
end
end
end
%constraint 3: We assume machines are initially in off state
for i = 1:M
W(i,1,1) == 1;
end
%constraint 4: Governs transition of machines, TA is set of allowed transition for the particular state. For the next period p+1, machine needs to follow this
for i = 1:M
for p = 1:P-1
for s = 1:S
for k = TA{s}
W(:,:,:) <= sum (W(:,k,p+1));
end
end
end
end
%constraint 5: Only one job can be processed at a time, sum over j for machine i in any period is 1
for i = 1:M
for p = 1:P
for j = 1:J
sum (X(:,j,:)) <= 1;
end
end
end
%constraint 6: All jobs are started and are only allowed to start once (WHY) sum over machines and periods, Y (machine starts processing) is 1 hence one of the machines in any period will be processing at least one and exactly one jobs
for j = 1:J
sum(sum(Y(i,:,p))) == 1;
end
%constraint 7: Does not allow pre-emption. When machine starts processing a job, it will take the proc. time given for the job
for i = 1:M
for j = 1:J
for p = 1:P - ProcTime(j)
for k = p:p + ProcTime(j)
Y(i,j,p) .* ProcTime(j) <= sum (X(:,:,k));
end
end
end
end
%constraint 8: Max peak energy will be atleast equal to sum of energy consumption over all the machines, over all the states, over the period k, divided by number of rolling window periods which is considered for determining the max demand during the day by utiliti company
for p = 1:P - RollingWindow
for k = p:p + RollingWindow
(1/RollingWindow) * (sum(sum(sum(W(i,s,k) .* EnergyDemand(s))))) <= d;
end
end
%constraint 9: For all states that require defined time period to finish, if machine was in same state in previous period as it is now, the transition did not occur
for i = 1:M
for s = TransitionStates
for p = 2:P
W(:,:,p) - W(:,:,p-1) <= Z(:,:,p);
end
end
end
%constraint 10: When machine enters a state which requires a specific time to complete, it will take that amount of time to complete which will be at most equal to sum of all such time periods elapsed by the particular machine, in that state
for i = 1:M
for s = TransitionStates
for p = 1: P - TransitionTime(s)
for t = p: p + TransitionTime(s)
TransitionTime(s) .* Z(i,s,p) <= sum (W(:,:,t))
end
end
end
end
%constraint 11: Turns off the machine by end of timespan by selecting an element in the set : Off, and puts machine in that state by specifying value of binary variable W as 1
for i = 1:M
for k = OffState
for p = P - TransitionTime(5)
sum (W(:,k,:)) == 1;
end
end
end
%constraint 12: Calculates earliness penalty
for j = 1:J
EarlyPenalty(:) >= DueDates(:) - CompletionTime(:);
end
%constraint 13: Calculates tardiness penalty
for j = 1:J
TardinessPenalty(:) >= CompletionTime(:) - DueDates(:);
end
%constraint 14: Calculates time it took to finish particular job by the machine, over the time period
for j = 1:J
CompletionTime(j) == sum (Y(:,:,p) .* p) + ProcTime(j);
end
cvx_end
%% debugs
%
% disp (d);
% for k = OffState
% for p = 1:P-TransitionTime{5}
% disp (p);
% disp (k);
% end
% end
%
% for p = P-TransitionTime{5}
% disp (p);
% end
%
%
% for s = TransitionStates
% disp(s);
% end
% for i =1:M
% for s = TransitionStates
% for p = 1: P - TransitionTime{s}
% for t = p: p + TransitionTime{s}
% disp (t);
% disp (p);
% disp (s);
% disp (i);
% end
% end
% end
% end
% for p = 1:P-1
% for s = 1:S
% for k = TA{s}
% disp (k);
% disp (s);
% disp (p);
% end
% end
% end
% for j = 1:J
% for p = 1:P-ProcTime(j)
% for k = p:p+ProcTime(j)
% disp (j);
% disp (p);
% disp (k);
% end
% end
% end
% for p = 1:P - RollingWindow
% for k = p:p + RollingWindow
% disp (p);
% disp (k);
% disp(RollingWindow);
% end
% end
% for p = P-TransitionTime(5)
% disp (p);
% end
% for s = TransitionStates
% for p = 2:P
% disp (s);
% disp (p);
% end
% end
}