ga optimization share variable between constraint and objective function to reduce computational time

I am working on an optimization problem using GA optimization. I have a variable involved in both the objective function adn the constraint function. This variable is days, and it is the constraint for the problem as well as it contributes to the cost, which is the objective function. The function that calculates the days variable is a very time consuming function, and having to calculate it twice it making the program take near double what it should. Is there any way for me to be able to only have this variable calculate once?

The calculation for days involves using the optimization variables within it, so I don't think I can just calculate it outside of the function since it changes as the optimization variables change.

Here is my code, there is quite a few side functions that I am not sharing here, but I can if absolutely necessary. I don't think the should entirely be required, the code works as is, it is just slower than I want, and I believe this will reduce that.

nvar = 2;
upperBounds = [10,10];
lowerBounds = [1,1];
intcon = 1:nvar;

addpath ../HelperFunctions/
addpath ../Welding/

rng default % for reproducibility
global inputs unitCostLaborFile unitCostFile weldSpeedsFile;

inputs = inputCaller();
[unitCostLaborFile, unitCostLaborPath] = uigetfile('*.xlsx','Select file with Welding Unit Costs for labor');
[unitCostFile, unitCostPath] = uigetfile('*.xlsx','Select file with Welding Unit Costs for equipment');
[weldSpeedsFile, weldSpeedsPath] = uigetfile('*.xlsx','Select file with Welding Speeds from WPS');

% Set nondefault solver options
options = optimoptions(@ga,'PlotFcn',{@gaplotbestf, @gaplotstopping});
options.PopulationSize = 10;

% Solve
fun = @objectiveFcn;
[solution,objectiveValue,exitflag,output,population,scores] = ga(fun,nvar,[],[],[],[],lowerBounds,...
    upperBounds,@constraintFcn,intcon,options);

% Clear variables
clearvars options

fill = solution(1);
cap = solution(2);
cost = objectiveValue;

disp(population(1:10,:));
disp(scores(1:10));
print(fill,cap,cost,inputs,weldSpeedsFile);

function f = objectiveFcn(optimInput)
% % Example:
% % Minimize Rosenbrock's function
% % f = 100*(y - x^2)^2 + (1 - x)^2
% 
% % Edit the lines below with your calculation
% x = optimInput(1);
% y = optimInput(2);
% f = 100*(y - x^2)^2 + (1 - x)^2;

global inputs unitCostLaborFile unitCostFile weldSpeedsFile;

% Input Variables
startDateInput = str2num(inputs{1});
endDateInput = str2num(inputs{2});
pipeSize = str2double(inputs(3));
pipeThickness = str2double(inputs(4));
welds = str2double(inputs(5));
fillPasses = str2double(inputs(6));
capPasses = str2double(inputs(7));
workHours = str2double(inputs(8));
travelTime = str2double(inputs(9));
weight = str2double(inputs(10));
overhang = str2double(inputs(11));
qaqc = str2double(inputs(12));
useTackRig = strcmp(inputs(13), 'true');
paidHours = ceil(workHours + travelTime);
startDate = datetime(startDateInput(1),startDateInput(2),startDateInput(3));
endDate = datetime(endDateInput(1),endDateInput(2),endDateInput(3));
formatIn = 'dd-mmm-yyyy';
deadline = datenum(endDate) - datenum(startDate);

% Cost Variables
welderCost = calculateDailyWage(unitCostLaborFile,"PL_WELDR",paidHours);
helperCost = calculateDailyWage(unitCostLaborFile,"PL_WH1",paidHours);
laborerSkilledCost = calculateDailyWage(unitCostLaborFile,"PL_LAB1",paidHours);
stabberCost = calculateDailyWage(unitCostLaborFile,"PL_ST",paidHours);
spacerO16Cost = calculateDailyWage(unitCostLaborFile,"PL_SO16",paidHours);
spacerU16Cost = calculateDailyWage(unitCostLaborFile,"PL_SU16",paidHours);
operatorCost = calculateDailyWage(unitCostLaborFile,"PL_OP1",paidHours);
foremanCost = calculateDailyWage(unitCostLaborFile,"PL_FM1",paidHours);
strawBossCost = calculateDailyWage(unitCostLaborFile,"PL_SB1",paidHours);
clampmanCost = calculateDailyWage(unitCostLaborFile,"PL_CL",paidHours);
qaqcCost = calculateDailyWage(unitCostLaborFile,"QA_SR1",paidHours);
pipelayerModel = determinePipelayer(weight, overhang);
equipmentMap = ExcelToMapObj(unitCostFile,1:3);
pipelayerCost = workHours * equipmentMap(pipelayerModel).UnitCost;
tackRigCost = workHours * equipmentMap("RTCKRIG").UnitCost;
truckPickupCost = paidHours * equipmentMap("PU02").UnitCost;
truckCrewCabDeckCost = paidHours * equipmentMap("CC03").UnitCost;
clampCost = workHours * equipmentMap("RCLMP20").UnitCost;
LOA = 175;

% Determine number of spacers/type required
if(pipeSize < 16)
    spacers = 1;
    spacerCost = spacerU16Cost;
else
    spacers = 2;
    spacerCost = spacerO16Cost;
end

% Determine the Number of Preheat Laborers Required
if(month(startDate) >= 10 || month(startDate) <= 4)
   preheatLaborers = 3; 
else
   preheatLaborers = 1;
end

if(useTackRig == true)
    tackRig = 1;
    pipelayers = 1;
else
    tackRig = 0;
    pipelayers = 2;
end

otherCrew = 4 + spacers + pipelayers + tackRig + qaqc; % straw,foreman,stabber,clampman = 4

hoursWelding = paidHours - 1.75; %1.75 = 0.5 lunch, 0.5 coffee, 0.75 setup/morning meetings
averageTimeBetweenPasses = 3; %min

% optimization variables 
fill = optimInput(1);
cap = optimInput(2);

welderPairs = fill + cap + 2; %assuming one root pass pair and one hot pass pair

% variables that rely on welders
helpers = welderPairs*2;
laborers = round(welderPairs*2/3 + preheatLaborers);
crew = welderPairs*4 + laborers + otherCrew;

% Determine the number of trucks needed for the crew (note welders and
% helpers bring their own trucks and the foreman,strawboss and qaqc take 
% pickup 4x4 trucks):P
crewNeedTruck = otherCrew - 2 - qaqc + laborers;
crewCabTrucks = ceil(crewNeedTruck/3);
pickupTrucks = 2 + qaqc; % foreman and strawboss

% the amount of days based on the amount of welders and how fast they are
% going 
days = (welderTimeCons(welds,fill,cap,fillPasses, capPasses,averageTimeBetweenPasses,weldSpeedsFile, pipeThickness, pipeSize)/hoursWelding);

% Cost function and inequality function
f = days*(welderPairs*2 * welderCost + helpers * helperCost + laborers * laborerSkilledCost...
    + foremanCost + strawBossCost + stabberCost + spacers*spacerCost + qaqc*qaqcCost...
    + clampmanCost + clampCost + pipelayers*pipelayerCost + tackRig*tackRigCost + (pipelayers+tackRig)*operatorCost...
    + crew*LOA + crewCabTrucks*truckCrewCabDeckCost + pickupTrucks*truckPickupCost);
end

function [c,ceq] = constraintFcn(optimInput)
% % Example:
% % Constrain a solution to the region
% % x^2 + y^2 <= 5
% % x^2 + y^2 >= 2
% % y = x^3
% 
% % Edit the lines below with your calculation
% % Note, if no inequality constraints, specify c = []
% % Note, if no equality constraints, specify ceq = []
% x = optimInput(1);
% y = optimInput(2);
% c(1) = x^2 + y^2 - 5;
% c(2) = 2 - x^2 - y^2;
% ceq = y - x^3;
fill = optimInput(1);
cap = optimInput(2);

global inputs unitCostLaborFile unitCostFile weldSpeedsFile;

% Input Variables
startDateInput = str2num(inputs{1});
endDateInput = str2num(inputs{2});
pipeSize = str2double(inputs(3));
pipeThickness = str2double(inputs(4));
welds = str2double(inputs(5));
fillPasses = str2double(inputs(6));
capPasses = str2double(inputs(7));
workHours = str2double(inputs(8));
travelTime = str2double(inputs(9));
paidHours = ceil(workHours + travelTime);
startDate = datetime(startDateInput(1),startDateInput(2),startDateInput(3));
endDate = datetime(endDateInput(1),endDateInput(2),endDateInput(3));
deadline = datenum(endDate) - datenum(startDate);

hoursWelding = paidHours - 1.75; %1.75 = 0.5 lunch, 0.5 coffee, 0.75 setup/morning meetings
averageTimeBetweenPasses = 3; %min

days = (welderTimeCons(welds,fill,cap,fillPasses,capPasses,averageTimeBetweenPasses,weldSpeedsFile, pipeThickness, pipeSize)/hoursWelding);

c = days - deadline;
ceq = [];
end

Thank you

Answer 1

[Walter Roberson]

You can use memoize() https://www.mathworks.com/help/matlab/ref/memoize.html

Note that it is a mistake to assume that the cost function is always going to be run before executing the nonlinear constraint function for the same inputs, and it is also a mistake to assume that the nonlinear constraint function is always going to be run before executing the cost function for the same inputs. ga() and fmincon() are permitted to run the cost function and the constraint function in any order they care to.

In particular, fmincon() first runs the cost function N+1 times for a function with N variables, without running the nonlinear constraints at all. It does this in order to do an initial gradient estimation. ga() does not do gradient estimation, but it still needs N+1 evaluations in order to create an initial simplex.

After the initial phase, at any point ga() and fmincon() might run the nonlinear constraint function a number of times in a row, trying to find the inside of the nonlinear inequality constraint boundary, and trying to find the edge of the nonlinear equality constraint boundary. Once it has decided it understands the boundary well enough, then the value it passes to the objective function will not necessarily be the same as the last value it passed to the nonlinear constraint function: the functions are permitted to look around the surface if they think they are in a safe region, only bothering to check back with the nonlinear constraint if the new point is one of the N+1 best nearby points and so potentially useful for updating the gradient estimation...

Because of these factors, calculating using a shared variable is riskier than it sounds at first, so memoization instead is safer.

[Jadon Latta]

Thank you so much, I wasn't aware of memoization, this is a very useful tool. I appreciate your help

Answer 2

[Sean de Wolski]

All of your constraints are simple arithmetic. I think you could solve this much more quickly and effectively as a mixed integer linear programming problem with intlinprog. Take a look at this example: Factory, Warehouse, Sales Allocation Model: Solver-Based - MATLAB & Simulink (mathworks.com)

[Jadon Latta]

Thanks Sean, this is actually only the first phase of the overall project, we know that there are other quicker and more effective methods to do something this simple, but we are intending to significantly complicate the equation over the next few months, as well as in the very long term implement machine learning into the software, hence the usage of GA