How to check if my linear constraints and bounds are infeasible & reduce execution time?

Hello,

I'm making a code where it creates an optimized engine usage schedule based on engine's reliability using genetic algorithm solver. The following is the code that I've made:

tic
clear
% Create arrays to store total Hour & daily result
totalHourWeekday = zeros(5,13);
result = zeros(5,13);

% Calculate the maximum engine by adding all of the engine's status
engineStatus = ones(1,13);
maxEngine = sum(engineStatus);

% Initialize optimization problem
rel = zeros(1,13);
decision = optimproblem("ObjectiveSense","max");
engine = optimvar('engine',13,'Type','integer','LowerBound',0,'UpperBound',1);

% Define days
for day = 1:5
    % calculate maximum duration
    maxDuration = 24*max(day);
    % Initialize demand
    for demand = 1:5
        disp(['Day: ',num2str(day),' Demand: ',num2str(demand)])
        switch demand
            case 1
                hourDemand = 3;
                engineDemand = 9;
                % Calculate reliability based on the sum between the hour
                % demand & previous hour
                if day == 1
                    timeCalc = hourDemand + totalHourWeekday(demand,:);
                else
                    timeCalc = hourDemand + result(day-1,:);
                end
                rel = exp(-(4.44*10^-5*timeCalc));
                relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));
               
                % Set reliability as objective function
                decision.Objective = relSys;
                
                % Initialize constraint #1: engines should be equal to the
                % demand
                engineMaintOrNot = sum(engine*engineStatus);
                decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;

                % Initialize constraint #2: engine's hour should be less
                % than max duration
                if day == 1
                    decision.Constraints.durationConstraint = engine*totalHourWeekday(demand,:) <= maxDuration;
                else
                    decision.Constraints.durationConstraint = engine*result(day-1,:) <= maxDuration;
                end
                % Solve optimization problem
                sol = solve(decision);
                
                % Add the hour demand to the chosen engines
                for i = 1:13
                    if sol.engine(i) == 1
                        totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
                    end
                end
    
                case 2
                    hourDemand = 2;
                    engineDemand = 6;
    
                    timeCalc = hourDemand + totalHourWeekday(demand-1,:);
                    rel = exp(-(4.44*10^-5*timeCalc));
                    relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));

                    decision.Objective = relSys;

                    engineMaintOrNot = sum(engine*engineStatus);
                    decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
                    
                    decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
                    sol = solve(decision);

                    for i = 1:13
                        if sol.engine(i) == 1
                            totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
                        end
                    end
                case 3
                    hourDemand = 5;
                    engineDemand = 10;

                    timeCalc = hourDemand + totalHourWeekday(demand-1,:);
                    rel = exp(-(4.44*10^-5*timeCalc));
                    relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));

                    decision.Objective = relSys;

                    engineMaintOrNot = sum(engine*engineStatus);
                    decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
                    
                    decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
                    sol = solve(decision);

                    for i = 1:13
                        if sol.engine(i) == 1
                            totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
                        end
                    end
                case 4
                    hourDemand = 8;
                    engineDemand = maxEngine-1;

                    timeCalc = hourDemand + totalHourWeekday(demand-1,:);
                    rel = exp(-(4.44*10^-5*timeCalc));
                    relSys = 1-((1-engine(1)*rel(1,1))*(1-engine(2)*rel(1,2))*(1-engine(3)*rel(1,3))*(1-engine(4)*rel(1,4))*(1-engine(5)*rel(1,5))*(1-engine(6)*rel(1,6))*(1-engine(7)*rel(1,7))*(1-engine(8)*rel(1,8))*(1-engine(9)*rel(1,9))*(1-engine(10)*rel(1,10))*(1-engine(11)*rel(1,11))*(1-engine(12)*rel(1,12))*(1-engine(13)*rel(1,13)));

                    decision.Objective = relSys;

                    engineMaintOrNot = sum(engine*engineStatus);
                    decision.Constraints.engineDemandConstraint = engineMaintOrNot == engineDemand;
                    
                    decision.Constraints.durationConstraint = engine*totalHourWeekday(demand-1,:) <= maxDuration;
                    sol = solve(decision);
    
                    for i = 1:13
                        if sol.engine(i) == 1
                            totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + hourDemand;
                        end
                    end
                case 5
                    % Special case
                    for i = 1:13
                        if engineStatus(1,i) == 1
                            totalHourWeekday(demand,i) = totalHourWeekday(demand,i) + 6;
                        end
                    end
        end
            % Add the total hour from all hour demands within the day 
            if day == 1
                result(day,:) = sum(totalHourWeekday);
            else
                result(day,:) = result(day-1,:) + sum(totalHourWeekday);
            end
    end
end

The code ran well until it reached 4th day. When calculating 4th day's 1st demand, it gave me the following error:

I guessed that, since it gave the infeasible linear constraints and bounds, an error was detected on line 57. That's why I am wondering how to solve this problem or how to check for linear constraints and bounds infeasibility?

Furthermore, is there any way to reduce the code's execution time?

Thank you in advance

[Torsten]

I'm surprised that you think you can solve your problem for each day separately and not "as one" for the complete week because of the couplings that arise.

[Muhammad Kahari Pranata]

@Torsten The original plan was to optimize which engines should be chosen for the hour demands based on each of the engine's reliability and I thought the only way to approach this was to "calculate" the chosen engines per demand per day. As an update, I've managed to solve the error by resetting the totalHourWeekday array every new day instead of outside the day for-loop. But, now that you have mentioned that, if you don't mind, could you tell me how to approach this optimization problem each week instead of each day?

[Torsten]

I don't know the background of your problem, but usually - if you make schedules for a certain time span - you cannot do this successively (day by day), but for the complete period because decisions taken at day i influence decisions on days before and after day i.

[Muhammad Kahari Pranata]

@Torsten I forgot to add some further information regarding this program. This code is just about one part of the whole project. The other part is to make sure that when one of the engines reached a certain operation hour, it must be maintained for certain period of days and make sure that maintenance happens only to one engine at a time (so multiple engines can't be maintained simultaneously). So the final output of this program is to generate a maintenance schedule that shows, for each engines, the day of the maintenance and the operating hour when the maintenance is executed. Can I still do it per week? or do I have no choice but to do it per day?

Answer 1

[Nipun]

Hi Kahari,

I understand that you are receiving an out of bounds error when indexing in a data structure in MATLAB and intend to know the root cause of this error.

You can create a filter to investigate the infeasibilities of the linear programming constraints in MATLAB. For more information on creating such filters, refer to the following MathWorks documentation: https://www.mathworks.com/help/optim/ug/investigate-linear-infeasibilities.html

Hope this helps.

Regards,

Nipun

[Muhammad Kahari Pranata]

Thank you for the answer, @Nipun. I've actually solved the problem by resetting the totalHourWeekday every new day instead of initializing it only once. But I will look into that if I encountered the error again. Now, the only problem is to reduce the code's execution time. Do you know how to reduce it?