Why is the following code working fine when executed with a for loop but showing an error when it comes to parfor?

The following code shows an "Index exceeds matrix dimensions error when run with parfor but runs okay when executed with just a for loop. mpc is a struct, in the code.

What am I doing wrong here and how do I fix it?

 %Beginning
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
mpc = loadcase('case9');
refbus = find(REF == mpc.bus(:,BUS_TYPE),1); 
nb = size(mpc.bus,BUS_I);
ng = size(mpc.gen,GEN_BUS);
nl = size(mpc.branch,F_BUS);
[Vx_idx,Vy_idx,Pf_idx,Qf_idx,Pt_idx,Qt_idx,Pg_idx,Qg_idx,Pf1_idx,Qf1_idx,Pt1_idx,Qt1_idx] = assign_idx(nb, nl, ng);
[gamma, theta, alpha1, alpha2, ~, termcrit] = flexa_constants; 
Aeq = [];
beq = [];
lb = [];
ub = [];
Ev_idx =0;
Evv_idx =0;
Ef_idx = 0;
Et_idx = 0;
optimval = 5216.028;
count = 1;
rel_val = 0;
relerror = 1;
tempx0 = 0;
cost_final = 0;
exitflag = 0;
output = 0;
ef = 0;
costs = zeros(nb,1);
Mtemp = 0;
Mvector = 0;
c = 0;
ceq = 0;
xsol = 0;
ef = zeros(nb,1);
X = cell(nb,1);
posf = 0;
post = 0;
pos = 0;
x_ii = 0;
%%Initial Guess for x
x_l = 2*nb+8*nl+4*ng;
Ev_idx = x_l+1:x_l+4*nb;
x_l = Ev_idx(end);
Evv_idx = x_l+1:x_l+nb;
x_l = Evv_idx(end);
Ef_idx = x_l+1:x_l+4*nl;
x_l = Ef_idx(end);
Et_idx = x_l+1:x_l+4*nl;
x_l = Et_idx(end);
x0 = zeros(x_l,1);
x0(Vx_idx) = 1;
x0(Pg_idx) = 1.05;
%%fmincon
% while relerror >= termcrit
while count == 1
    tempx0 = x0;
    H = sparse(x_l,x_l);    %Quadratic term coefficient
    h = sparse(x_l,1);  %Linear term coefficient
    g1_idx = Pg_idx(1:3);
    g2_idx = Pg_idx(4:6);    
    q1_idx = Qg_idx(1:3);
    q2_idx = Qg_idx(4:6);    
    A0 = [-eye(nb) zeros(nb, x_l-nb)];
    b0 = zeros(nb,1);
    A1 = zeros(2,x_l);
    A1(1,g1_idx) = -1;
    A1(2,g2_idx) = -1;
    A = [A0;A1];
    b1 = [-sum(mpc.bus(:,PD))/100;-sum(mpc.bus(:,PD))/100];
    b = [b0;b1];    
    %
    Vlb = zeros(x_l,1);
    Vub = zeros(x_l,1);
    Vlb(1:x_l) = -Inf;
    Vub(1:x_l) = Inf;
    Vlb(Vx_idx) = 0.9;
    Vub(Vx_idx) = 1.1;
    for gg = 1:ng
        Vlb(Pg_idx(gg)) = mpc.gen(gg,PMIN)/100;
        Vlb(Pg_idx(ng+gg)) = Vlb(Pg_idx(gg));
        Vub(Pg_idx(gg)) = mpc.gen(gg,PMAX)/100;
        Vub(Pg_idx(ng+gg)) = Vub(Pg_idx(gg));
        Vlb(Qg_idx(gg)) = mpc.gen(gg,QMIN)/100;
        Vlb(Qg_idx(ng+gg)) = Vlb(Qg_idx(gg));
        Vub(Qg_idx(gg)) = mpc.gen(gg,QMAX)/100;
        Vub(Qg_idx(ng+gg)) = Vub(Qg_idx(gg));
    end
    lb = Vlb;
    ub = Vub;    
    const = 0; 
    for i = 1:ng
        H(Pg_idx(i),Pg_idx(i)) = mpc.gencost(i,COST);
        h(Pg_idx(i),1) = mpc.gencost(i,COST+1);
        const(i) = mpc.gencost(i,COST+2);
    end    
    gradF = 2*H*100*x0 + h*100;
    mag = norm(gradF);
    epsilon = gamma*alpha1*min(alpha2,1/mag);    
    Aslack = zeros(1,x_l);
    Aslack(1,Vy_idx(refbus)) = 1;
    Aeq = [Aslack]; beq = 0;
    % epsilon = 1e-3;    
    fun = @(x)(100*x)'*H*(100*x) + (100*x)'*h + sum(const)...
        + power(norm(x(g1_idx)-x(g2_idx)),2)...
        + power(norm(x(g1_idx)-x(g2_idx)),2)...
        + power(norm(x(Pf_idx)-x(Pf1_idx)),2) + power(norm(x(Pt_idx)-x(Pt1_idx)),2)...
        + power(norm(x(Qf_idx)-x(Qf1_idx)),2) + power(norm(x(Qt_idx)-x(Qt1_idx)),2)...
        + power(norm(x(Ev_idx)),2) + power(norm(x(Evv_idx)),2) + power(norm(x(Ef_idx)),2) + power(norm(x(Et_idx)),2);
    options = optimoptions(@fmincon,'Algorithm','interior-point','ConstraintTolerance', 1e-4,'MaxIterations',1e6,'MaxFunctionEvaluations',1e6,'StepTolerance',epsilon);
    g = zeros(nb,1);
    g(mpc.gen(:,GEN_BUS)) = mpc.gen(:,GEN_BUS);
    parfor bus_idx = 1:nb
        [x{bus_idx},fval(bus_idx),exitflag(bus_idx),output(bus_idx),lambda,grad,hessian] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,@(x)nonlconstraints_editpar1(x,x0,Ev_idx,Evv_idx,Ef_idx,Et_idx,bus_idx),options);
        ef(bus_idx) = exitflag(bus_idx);
        costs(bus_idx) = fval(bus_idx);           
        posf = find(bus_idx==mpc.branch(:,1));
        post = find(bus_idx==mpc.branch(:,2));
        pos = union(posf,post);        
        if g(bus_idx)~=0
            x_ii = [Vx_idx(bus_idx);Vy_idx(bus_idx);Pf_idx(pos);...
                Qf_idx(pos);Pt_idx(pos);Qt_idx(pos);Pg_idx(bus_idx);...
                Qg_idx(bus_idx);Pf1_idx(pos);Qf1_idx(pos);Pt1_idx(pos);...
                Qt1_idx(pos)];
        else
            x_ii = [Vx_idx(bus_idx);Vy_idx(bus_idx);Pf_idx(pos);...
                Qf_idx(pos);Pt_idx(pos);Qt_idx(pos);...
                Pf1_idx(pos);Qf1_idx(pos);Pt1_idx(pos);...
                Qt1_idx(pos)];
        end
        xsol = x0(x_ii) + gamma*(x{bus_idx}(x_ii) - x0(x_ii));
        X{bus_idx} = xsol;
    end
    count = count + 1;
end

Error:
>> optim_par_nogrp
Error using optim_par_nogrp>(parfor supply)
Index exceeds matrix dimensions.
Error in optim_par_nogrp (line 137)
  parfor bus_idx = 1:nb

[Jan]

Just a remark: The power(x,2) operator is much more expensive than x.^2.

Is this a script of a function?

[Viswanath Hariharan]

Yes it is.

[ADragon]

Hi Viswanath, I would take everything inside of the parfor loop and create a function. Then call that function using the indexer.

parfor bus_idx = 1:nb
    X{bus_idx} = myfunction(bus_idx,...);
end

Do you get the same behavior? I have found that the parfor loop does really well to populate a single array given the index as the looping variable. Since you have multiple arrays inside parfor, it may mess up how the parfor function sets up the different threads. If you isolate the array you want to parallel compute for (as the above code) it may help solve this problem.