Using parallel computing for problem-based optimization

Hi there,

I want to speed up my program which uses problem-based optimization by parallel computing, I set 'UseParallel' to be true as optimization options. However, it came up with following error.

RCfun3R2C (line 4)
Index in position 2 exceeds array bounds.

Error in generatedObjective (line 30)

Error in optim.problemdef.OptimizationProblem/compileObjectives>@(x)objhandle.Value(x,extraParams.Value)

Error in finDiffEvalAndChkErr

Error in parfinitedifferences

Error in parfinitedifferences

Error in computeFinDiffGradAndJac

Error in sfdnls (line 54)
computeFinDiffGradAndJac(x,funfcn,confcn,valx, ...
    
Error in snls (line 178)
[A,findiffevals] = sfdnls(xcurr,fvec,Jstr,group,alpha,funfcn,l,u,...
    
Error in lsqncommon (line 164)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...
    
Error in lsqnonlin (line 262)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...
    
Error in optim.problemdef.OptimizationProblem/callSolver

Error in optim.internal.problemdef.ProblemImpl/solveImpl

Error in optim.problemdef.OptimizationProblem/solve

Error in RC_calibration (line 42)
[rsol,~] = solve(prob,r0, 'Options', options);

Error in Main (line 17)
RC_parameters = RC_calibration(data_train, data, need_plot);

My objective is RCfun3R2C, and I want to solve optimization problem to find the best model parameters:

function T = RCfun3R2C(r)
% weather parameters
parameter_RCfun = theGlobalData().parameter_RCfun;
temp_381A = parameter_RCfun(:, 1);
temp_amb = parameter_RCfun(:, 2);
temp_wall_out = parameter_RCfun(:, 3);
pHVAC_381A = parameter_RCfun(:, 4);
Q_sol = parameter_RCfun(:, 5);
Q_int = parameter_RCfun(:, 6);


% Building parameters
R1 = r(1);
R2 = r(1);
Rwin = r(2);
Cz = r(3);
Cw = r(4);
ratio_sol = r(5);
% Wall_absorb = r(7);
% Win_trans = r(8);

% Build 5R3C model
A = [ -(R1+R2)/(Cw*R1*R2), 1/(Cw*R1); ... 
    1/(Cz*R1), -(R1+Rwin)/(Cz*R1*Rwin)];
B = [0; 1/Cz]; 
E = [1/(Cw*R2), 1/Cw, 0; ...
    1/(Cz*Rwin), 1/Cz, 1/Cz];
C = [0, 1];
D = 0;

% Discretize
dt = 60*5;
Ad = expm(A*dt);
Bd = inv(A)*(Ad-eye(2))*B;
Ed = inv(A)*(Ad-eye(2))*E;

% Solar radiation
Q_sol = Q_sol*ratio_sol;

% Closed loop
x0 = [temp_wall_out(1), temp_381A(1)];
x = zeros(2, length(temp_381A));
x(:, 1) = x0;
u = pHVAC_381A';
w = [temp_amb'; Q_sol'; Q_int'];
for t = 1:length(temp_381A)-1  
    x(:,t+1) = Ad*x(:,t) + Bd*u(t) + Ed*w(:,t);
end
T = x(2,:);
end

The code for optimization problem setup is showing below:

function RC = RC_calibration(Training_data, Testing_data, need_plot)
%% Data process
Training = Training_data;
Testing = Testing_data;

%% Calibration
% Global variable
parameter_training = [Training.temp_381A,  Training.temp_amb, ...
    Training.temp_wall_out, Training.pHVAC_381A, ...
    Training.Q_sol, Training.Q_int];
theGlobalData('parameter_RCfun', parameter_training);
theGlobalData('temp_381A', parameter_training(:, 1));

% Calibrate model
yvalstrue = parameter_training(:, 1)';
r = (1:5);
T = RCfun3R2C(r);

% Optimization problem
r = optimvar('r',5,"LowerBound",[0.2,...  % R1 = Exterior wall thermal resistance
    0.01,... % Window thermal resistance
    5e6,... % Zone heat capacity
    1.5e7,... % Exterior wall heat capacity
    0],... % Solar radiation modification coefficient
    "UpperBound",[0.6,...
    0.03,...
    3e7,...
    1e8,...
    2]);
myfcn = fcn2optimexpr(@RCfun3R2C,r);
obj = sum((myfcn - yvalstrue).^2);
prob = optimproblem("Objective",obj);
show(prob)

% Solve Problem
r0.r = [0.4,...
    0.02,...
    1e7,...
    3e7,...
    1/30];
options = optimoptions(prob,'UseParallel',true);
[rsol,~] = solve(prob,r0, 'Options', options);
disp(rsol.r)
temp_calibration = RCfun3R2C(rsol.r);
MAPE_calibration = mean(abs((temp_calibration-yvalstrue)./yvalstrue));

% Validation
parameter_testing = [Testing.temp_381A, Testing.temp_amb, ...
    Testing.temp_wall_out, Testing.pHVAC_381A, ...
    Testing.Q_sol, Testing.Q_int];
theGlobalData('parameter_RCfun', parameter_testing);
temp_validation = RCfun3R2C(rsol.r);
MAPE_validation = mean(abs((temp_validation-Testing.temp_381A')./Testing.temp_381A'));

%% Output RC Parameters
RC = rsol.r;

%% Plot calibration and validatation results
if need_plot == 1
    figure
    plot(Training.time, temp_calibration, 'LineWidth', 1)
    hold on
    plot(Training.time, Training.temp_381A, 'LineWidth', 1)
    % plot(Training.time, Training.temp_381B, 'LineWidth', 2)
    plot(Training.time, Training.temp_amb, 'LineWidth', 1)
    set(gca, 'Fontsize', 14);
    set(gca,'fontsize', 14, 'ygrid','on','xgrid','on');
    datetick('x','mmm-dd HH:MM','keepticks');
    ylabel('temperature (°C)');
    % ylim([min(temp_amb),  max(max(temp_381A), max(temp_381B))]);
    title(['3R2C Training MAPE: ', sprintf('%g', round(MAPE_calibration*100, 2)), '%'])
    legend('Prediction', 'Measurement(381A)', 'Ambient temperature','FontSize', 14)

    figure
    plot(Testing.time, temp_validation, 'LineWidth', 1)
    hold on
    plot(Testing.time, Testing.temp_381A, 'LineWidth', 1)
    % plot(Testing.time, Testing.temp_381B, 'LineWidth', 2)
    plot(Testing.time, Testing.temp_amb, 'LineWidth', 1)
    set(gca, 'Fontsize', 14);
    set(gca,'fontsize', 14, 'ygrid','on','xgrid','on');
    datetick('x','mmm-dd HH:MM','keepticks');
    ylabel('temperature (°C)');
    % ylim([mxin(temp_amb),  max(max(temp_381A), max(temp_381B))]);
    title(['3R2C Testing MAPE: ', sprintf('%g', round(MAPE_validation*100, 2)), '%'])
    legend('Prediction', 'Measurement(381A)', 'Ambient temperature','FontSize', 14)
end
end

I searched lots of parallel computing issues but didn't find one relating with mine.

Could anyone please help me with it? I would really appreciate!

Best regards,

Wang

[Mario Malic]

Try

dbstop if error

RCfun3R2C (line 4)

Index in position 2 exceeds array bounds.

Inspect the array parameter_RCfun, calculated values might be wrong/NaN. Also, why do you use parentheses to index into structure, not sure if that's neccessarry.

function T = RCfun3R2C(r)
% weather parameters
parameter_RCfun = theGlobalData().parameter_RCfun;
temp_381A = parameter_RCfun(:, 1);
temp_amb = parameter_RCfun(:, 2);
temp_wall_out = parameter_RCfun(:, 3);

[Xuezheng Wang]

Hi Mario,

Thank you for the suggestions. I've tried but same error still happened.

I have checked the parameter_RCfun for several times, all calculations are right and there is no NaN.

For the theGlobalData(), it is used as a surrogate of global variable.

function G = theGlobalData(Name, Value)
persistent G_
if isempty(G_)
    G_.parameter_RCfun = [];
end
if nargin > 0
    G_.(Name) = Value;
end
G = G_;
end

What weird is that when I turn off parallel computing, the program could run with no error.

Answer 1

[Max Heimann]

If the program runs fine without parallel computing and errors out with parallel computing the error might be a dynamic array size. Usually matlab is fine with arrays changing size in e.g. loops. However when using parallel computing this is not allowed for various reasons. Tihs could lead to your index out of bonds error. I would suggest you check your code whether there are any dynamicly sized variables. (e.g. appending entries to a vector, adding an i-th element to a vector in a loop over i...).

The solution might be to preallocate those variables.

[Xuezheng Wang]

Hi Max! Thank you for the suggestions. I found the problem is related with my global variable function. It is solved by nested function.