Optimize matrices operations

Question

I'm optimizing via Genetic Algorithm a machine simulator (of a MultiheadWeigher machine) in order to solve the well known "Setup problem". I based all the code on matricies but with the multiproduct case I think there is still some inefficiency...

Code

function [f] = MHW(position_matrix, HA, HB, HC)

global W_best_tot
global Function
global n_b_global

nComb = size(position_matrix,1);
dim = size(position_matrix,2);

WL = 241;

%SIMULATIONS VARIABLES
alpha = 0.123;
nSim = 3000;
CellUncertainty = 0.5 * 10^(-6);

%// Define Product Hopper allocation
WT_A = 130; WL_A = 129;
WT_B = 30; WL_B = 29;
%HC = 2; 
WT_C = 90; WL_C = 89;

%// COMBINATION MATRIX

CombinantionMatrix_A = combn([0 1], HA);
CombinantionMatrix_B = combn([0 1], HB);
CombinantionMatrix_C = combn([0 1], HC);

if HA == 1
    CombinantionMatrix_A = CombinantionMatrix_A.';
end
if HB == 1
    CombinantionMatrix_B = CombinantionMatrix_B.';
end
if HC == 1
    CombinantionMatrix_C = CombinantionMatrix_C.';
end

CombinantionMatrix_A_Transp = CombinantionMatrix_A.';
CombinantionMatrix_B_Transp = CombinantionMatrix_B.';
CombinantionMatrix_C_Transp = CombinantionMatrix_C.';

% OBJECTIVE FUNCTION COST COEFFICIENTS
Cu_A = 0.03; c_p = 0.6; c_f = 734; c_l = 3.2;
Cu_B = 0.09;
Cu_C = 0.04;

[HopperWeight UncertainWeight] = Weight(nComb, dim, alpha, ...
    position_matrix, CellUncertainty);

HopperWeight_A = HopperWeight(:,1:HA);
HopperWeight_B = HopperWeight(:,HA+1:HA+HB);
HopperWeight_C = HopperWeight(:,HA+HB+1:HA+HB+HC);

UncertainWeight_A = UncertainWeight(:,1:HA);
UncertainWeight_B = UncertainWeight(:,HA+1:HA+HB);
UncertainWeight_C = UncertainWeight(:,HA+HB+1:HA+HB+HC);

W_best_tot = zeros(nComb, nSim);
W_best_ABC = zeros(nComb, 3, nSim);

for ii=1:nSim

    W_A = HopperWeight_A * CombinantionMatrix_A_Transp;
    W_B = HopperWeight_B * CombinantionMatrix_B_Transp;
    W_C = HopperWeight_C * CombinantionMatrix_C_Transp;

    W_Unc_A = UncertainWeight_A * CombinantionMatrix_A_Transp;
    W_Unc_B = UncertainWeight_B * CombinantionMatrix_B_Transp;
    W_Unc_C = UncertainWeight_C * CombinantionMatrix_C_Transp;

    [~,comb_A] = min(abs(W_Unc_A - WT_A),[],2);
    [~,comb_B] = min(abs(W_Unc_B - WT_B),[],2);
    [~,comb_C] = min(abs(W_Unc_C - WT_C),[],2);

    W_best_A = W_A(sub2ind_mod(size(W_A),1:size(W_A,1),comb_A.')).';
    W_best_B = W_B(sub2ind_mod(size(W_B),1:size(W_B,1),comb_B.')).';
    W_best_C = W_C(sub2ind_mod(size(W_C),1:size(W_C,1),comb_C.')).';

    W_best_ABC(:,:,ii) = [W_best_A W_best_B W_best_C];
    W_best_tot(:,ii) = sum(W_best_ABC(:,:,ii),2);

    [HopperWeight_2_A UncertainWeight_2_A] = Weight(nComb, HA, alpha, ...
        position_matrix(:,1:HA), CellUncertainty);
    [HopperWeight_2_B UncertainWeight_2_B] = Weight(nComb, HB, alpha, ...
        position_matrix(:,HA+1:HA+HB), CellUncertainty);
    [HopperWeight_2_C UncertainWeight_2_C] = Weight(nComb, HC, alpha, ...
        position_matrix(:,HA+HB+1:HA+HB+HC), CellUncertainty);

    idx = CombinantionMatrix_A(comb_A,:)~=0;
    HopperWeight_A(idx) = HopperWeight_2_A(idx);
    UncertainWeight_A(idx) = UncertainWeight_2_A(idx);
    idx = CombinantionMatrix_B(comb_B,:)~=0;
    HopperWeight_B(idx) = HopperWeight_2_B(idx);
    UncertainWeight_B(idx) = UncertainWeight_2_B(idx);
    idx = CombinantionMatrix_C(comb_C,:)~=0;
    HopperWeight_C(idx) = HopperWeight_2_C(idx);
    UncertainWeight_C(idx) = UncertainWeight_2_C(idx);

    clear HopperWeight_2_A HopperWeight_2_B HopperWeight_2_C ...
        UncertainWeight_2_A UncertainWeight_2_B UncertainWeight_2_C;

end


n_b = logical(W_best_tot >= WL);
n_bA = logical(reshape(W_best_ABC(:,1,:), nComb, nSim) >= WL_A);
n_bB = logical(reshape(W_best_ABC(:,2,:), nComb, nSim) >= WL_B);
n_bC = logical(reshape(W_best_ABC(:,3,:), nComb, nSim) >= WL_C);

n_b_global = sum(n_b .* n_bA .* n_bB .* n_bC, 2);

GAvg_A  = sum(reshape(W_best_ABC(:,1,:), nComb, nSim).*(reshape(...
            W_best_ABC(:,1,:), nComb, nSim) > WT_A),2)./sum(reshape(...
            W_best_ABC(:,1,:), nComb, nSim) > WT_A,2);
GAvg_B  = sum(reshape(W_best_ABC(:,2,:), nComb, nSim).*(reshape(...
            W_best_ABC(:,2,:), nComb, nSim) > WT_B),2)./sum(reshape(...
            W_best_ABC(:,2,:), nComb, nSim) > WT_B,2);
GAvg_C  = sum(reshape(W_best_ABC(:,3,:), nComb, nSim).*(reshape(...
            W_best_ABC(:,3,:), nComb, nSim) > WT_C),2)./sum(reshape(...
            W_best_ABC(:,3,:), nComb, nSim) > WT_C,2);

f = Cu_A .* GAvg_A + Cu_B .* GAvg_B + Cu_C .* GAvg_C + (nSim./...
        n_b_global) .* c_p + (c_f./n_b_global) + ((nSim - n_b_global)./...
        n_b_global) .* c_l;

Function = f;

Profiler

This is the Main Profiler and these are: the MHW function and the Weight function. The main problems are in the Weight since it's called 3000times for each GA generation...Above the Weight function code:

Weight Function Code

function [ HopperWeight UncertainWeight ] = Weight( nComb, H, alpha, AllComb, CellUncertainty )

HopperWeight = randn(nComb,H) * alpha;
HopperWeight = (1+HopperWeight) .* AllComb;
idx = HopperWeight<0;
random = randn(sum(idx(:)), 1);
HopperWeight(idx) = random .* ((1+alpha) .* AllComb(idx));

%ADD ERROR
RandomizeUncertainty = randn(nComb,H) * CellUncertainty;
UncertainWeight = abs((1+RandomizeUncertainty) .* HopperWeight);
end

Is there something that I am missing in order to optimize better the Weight function and in general the time consuming part in the MHW function? (If you need the GA launcher in order to try the MHW function use THIS code)

TIA