Infeasible Problem for gurobi 7.52 but incorrectly solved by 8.11 (likely formulation error)

Hello

When trying to solve my job shop scheduling problem, Gurobi 8.11 says its solved [Status: Solved
Optimal value (cvx_optval): +0] but optimal value is 0 while cvx supplied Gurobi 7.52 [Status: Infeasible
Optimal value (cvx_optval): +Inf] says the model is infeasible. I have formulated the exact same problem in AMPL and solved using CPLEX. Likely there is something wrong with the formulation, but I went through each and every constraint and loop that I could to check what is going on but could not find any errors (likely my MATLAB understanding has failed me) using my debugging techniques. I tried using other bundled solvers but they failed as well. I am attaching a reproducible codes without any comments but I can provide comments and expansions (or what the constraints should actually look like) in case that is needed.
{
P = 32;
S = 5;
M = 2;
J = 4;
DemandCharge = [10];
RollingWindow = [1];
ProcTime = [10 9 10 9];
TA{1} = [1 2];
TA{2} = [2 3 4];
TA{3} = [3 4 5];
TA{4} = [3 4 5];
TA{5} = [1 5];
ProcState = [4];
TransitionStates = [2 5];
OffState = [1 5];
TransitionTime = [0 5 0 0 1];
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];
EnergyDemand = [0 5 2 6 1];
cvx_begin
variable W(M,S,P) binary;
variable X(M,J,P) binary;
variable Y(M,J,P) binary ;
variable Z(M,S,P) binary;
variable d Nonnegative;
variable DemandCost Nonnegative;
variable MainCost Nonnegative;

    minimize ((DemandCost) + (MainCost);
    subject to
    
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

    for i = 1:M
        for s = ProcState
            for p = 1:P
            sum(X(:,j,:)) == sum(W(:,ProcState,:));
            end
        end
    end

    for i = 1:M
        for p = 1:P
            for s = 1:S
                sum (W(:,s,:)) == 1;
            end
        end
    end

    for i = 1:M
        W(i,1,1) == 1;
    end

    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

    for i = 1:M
        for p = 1:P
            for j = 1:J
            sum (X(:,j,:)) <= 1;
            end
        end
    end

    for j = 1:J
        sum(sum(Y(i,:,p))) == 1; 
    end

    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

    for p = 1:P - RollingWindow 
        for k = p:p + RollingWindow 
        (1/RollingWindow) * (sum(sum(sum(W(i,s,k) .* EnergyDemand(s))))) <= d;
        end
    end

    for i = 1:M
        for s = TransitionStates
            for p = 2:P
                W(:,:,p) - W(:,:,p-1) <= Z(:,:,p);
            end
        end
    end

    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

    for i = 1:M
        for k = OffState
            for p = P - TransitionTime(5)
            sum (W(:,k,:)) == 1; 
            end
        end
    end

    cvx_end
    }

Is the problem actually infeasible this way or there is problem with my solver or setup? How can I reformulate it? I have also attached picture of actual mathematical representation of the model.

I would greatly appreciate any input or feedback and even general advice.

Thank you in advance!!

1|393x500

You are missing a parenthesis in the minimize statement. Anyhow, I can’t execute youur program because the j used in sum(X(:,j,:)) == sum(W(:,ProcState,:)); is not defined, so it takes the value of the imaginary unit, which causes an error And the triple for loop it is in doesn’t seem to make any sense.

So more care is needed on your part.

Note that there have been dcoumented cases of Gurobi producing incorrect solutions under CVX – perhaps CVX/Gurobi interface bugs.

Sorry about that. I must have removed it while trying to remove my comments and some extra stuff to make it more clear. When I copy paste the code from here it gives me error as well but when I run the original which has exact same for loops and declarations, it runs without an error.

This is how the loop should be evaluated:

for i =1, s = 4 (which is defined as the value of ProcState), and p = 1,

sum over all values of j (from 1 to 4) on LHS

X(1,1,1) + X(1,2,1) + X(1,3,1) + X(1,4,1) + X(2,1,1) + X(2,2,1) + X(2,3,1) + X(2,4,1) == W(1,4,1) + W (2,4,1)

Since X and W are binary variables only two terms are allowed to be 1 on LHS at most.

This is how I was taught to write it for CVX to obtain this:

for i = 1:M
for s = ProcState
for p = 1:P
sum(X(:,j,:)) == sum(W(:,ProcState,:));
end
end
end

Hope it makes sense now

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
}

Mosek reports it as (primal) infeasible (when interpreted relative to your original problem). You have a big complicated model, which I don;'t know (haven’t checked) was entered (by you) correctly.

Just for the heck of it, I relaxed all the binary variables to be continuous 0 <= variable <= 1. and the problem is still reported infeasible.

I suggest you follow the advice in https://yalmip.github.io/debugginginfeasible/ .

Thanks again mark. Yes it is a big complicated model and it only gets bigger and more complicated. I have tried to make sure it was entered correctly by me but if you can maybe pick up just one constraint and confirm that I did enter it correctly, I would be very grateful. I will read the yalmpi article, thank you.