How can I write a matlab code for three index summation? Please help me.

16 views (last 30 days)
Where Q is a variable, B & C are parameters and Beta is a binary variable, whose value can be zero or one.
  4 Comments
Walter Roberson
Walter Roberson on 23 Feb 2017
I think you are saying that B & C are given by the user, but that Q is something that you have computed somehow?
For the purposes of the summation, it does not matter where the values came from.
Let's see:
Q must be 4 by 2 by 3
B must be 4 by 2
C must be 4 by 2 by 3
beta must be 4 by 2
You mention Beta as being a binary variable, but you use * before it, which would tend to imply it is somehow different than the multiplication of Q B C. Are we to understand that Q_dcp is a matrix to be algebraically multiplied by the matrix B_dc ? Or was that just an irrelevant oddity of notation?
John
John on 25 Feb 2017
Yes,these all(Q,B,C,Beta) are matrix and * is just used to separate binary variable from other values. Q is the value that we have to find.But in the start we will also put its value. I explain you by an example.Consider it is supply chain problem and d is showing number of plants,c=number of distributors and p=number of products.Now Q is the Quantity of different products that will move from each plant to different distributors.And it will be in form of matrix like d1 to c1,c2 or d2 to c1,c2. Index is basically explaining the Quantity of products move from each plant to distributors i.e. Q_dcp means Quantity of product p move from plants d to distributors c.And Beta is the binary variable matrix that will show the which plants or distributors are connected.

Sign in to comment.

Answers (4)

Jan
Jan on 23 Feb 2017
s = 0;
for d = 1:4
for c = 1:2
for p = 1:3
s = s + Q(d,c,p) * B(d, c) * C(d, c, p) * Beta(d, c);
end
end
end
  4 Comments
Jan
Jan on 23 Feb 2017
Edited: Jan on 23 Feb 2017
@John: I cannot guess, what Q, B and C is. Please do not state only, that it is not working, but explain the occurring problems. How could I find out, what is wrong with the answer and what you expect? See this code snippet as a suggestion, how such questions are solved usually and adjust it to your needs. Or provide enough information wuch that the wanted result can be reproduced.
John
John on 25 Feb 2017
These all(Q,B,C,Beta) are matrix and * is just used to separate binary variable from other values. Q is the value that we have to find.But in the start we will also put its value. I explain you by an example.Consider it is supply chain problem and d is showing number of plants,c=number of distributors and p=number of products.Now Q is the Quantity of different products that will move from each plant to different distributors.And it will be in form of matrix like d1 to c1,c2 or d2 to c1,c2. Index is basically explaining the Quantity of products move from each plant to distributors i.e. Q_dcp means Quantity of product p move from plants d to distributors c.And Beta is the binary variable matrix that will show the which plants or distributors are connected.Check attached file.

Sign in to comment.


Walter Roberson
Walter Roberson on 23 Feb 2017
If matrices are not involved, would this not be the same as
sum d from 1 to 4 of
sum c from 1 to 2 of
B(d,c) times beta(d,c) times
(sum p from 1 to 3 of
Q(d,c,p) * C(d,c,p) )
If so then
s = 0;
for d = 1 : 4
for c = 1 : 4
s = s + B(d,c) * beta(d,c) * sum(Q(d,c,:) .* C(d,c, :));
end
end
which then has the further potential for vectorization?
B_beta = B .* beta;
s = 0;
for d = 1 : 4
for c = 1 : 4
s = s + B_beta(d, c) * sum(Q(d,c,:) .* C(d,c, :));
end
end
and then to
B_beta = B .* beta;
sQC = sum(Q .* C, 3);
s = 0;
for d = 1 : 4
for c = 1 : 4
s = s + B_beta(d, c) * sQC(d,c);
end
end
which would lead to
B_beta = B .* beta;
sQC = sum(Q .* C, 3);
s = sum( sum( B_beta .* sQC ) );
or
B_beta = B .* beta;
sQC = sum(Q .* C, 3);
s = sum( B_beta(:) .* sQC(:) );
  3 Comments
John
John on 26 Feb 2017
Please do it with numeric values. Check the attached file for numeric values.
Thanks.
Walter Roberson
Walter Roberson on 26 Feb 2017
Using plain looping to establish a baseline of what the result should be:
%common initialization
Q_dcp = [200 300;250 0;0 200]; %from table
Beta_dcp = [1 1;1 0; 0 1]; %from table
B_dc = [2 1;3 2;1 4]; %from table
C_dcp = [5000 2000;6000 2000;1000 4500]; %from table
max_d = 3; %from table
max_c = 2; %from table
max_p = 1; %from table
%do the calculation
Q = Q_dcp;
Beta = Beta_dcp;
B = B_dc;
C = C_dcp;
s = 0;
for d = 1 : max_d
for c = 1 : max_c
for p = 1 : max_p
s = s + Q(d,c,p) .* B(d,c) .* C(d,c,p) .* Beta(d,c);
end
end
end
The result is s = 10700000
Using my more compact expression:
%common initialization
Q_dcp = [200 300;250 0;0 200]; %from table
Beta_dcp = [1 1;1 0; 0 1]; %from table
B_dc = [2 1;3 2;1 4]; %from table
C_dcp = [5000 2000;6000 2000;1000 4500]; %from table
max_d = 3; %from table
max_c = 2; %from table
max_p = 1; %from table
%do the calculation
Q = Q_dcp;
Beta = Beta_dcp;
B = B_dc;
C = C_dcp;
B_Beta = B .* Beta;
sQC = sum(Q .* C, 3);
s = sum( B_Beta(:) .* sQC(:) );
The result is s = 10700000, the same as the loop. This confirms that my compact code for the calculation does the same thing as the looping.
Using John BG's code:
%common initialization
Q_dcp = [200 300;250 0;0 200]; %from table
Beta_dcp = [1 1;1 0; 0 1]; %from table
B_dc = [2 1;3 2;1 4]; %from table
C_dcp = [5000 2000;6000 2000;1000 4500]; %from table
max_d = 3; %from table
max_c = 2; %from table
max_p = 1; %from table
%do the calculation
Q = Q_dcp;
Beta = Beta_dcp;
B = B_dc;
C = C_dcp;
S = zeros(max_d, max_d);
for k=1:1:p
S=S+Q(:,:,k).*B.*C(:,:,k)*Beta';
end
This gives the array result,
2600000 2000000 600000
4500000 4500000 0
3600000 0 3600000
If we then suggest that perhaps a couple of levels of summation was overlooked, and take sum(S(:)) then we arrive at
21400000
which is exactly twice the value arrived at by plain for loops.

Sign in to comment.


John BG
John BG on 25 Feb 2017
Edited: John BG on 25 Feb 2017
John
d=4;c=2;p=3;
Q=randi([-10 10],d,c,p)
B=randi([-10 10],d,c)
C=randi([-10 10],d,c,p)
beta=randi([-10 10],d,c)
S=zeros(d,d);
for k=1:1:p
S=S+Q(:,:,k).*B.*C(:,:,k)*beta';
end
if you find this answer useful would you please be so kind to mark my answer as Accepted Answer?
To any other reader, please if you find this answer of any help solving your question,
please click on the thumbs-up vote link,
thanks in advance
John BG
  1 Comment
Walter Roberson
Walter Roberson on 26 Feb 2017
John, I think you overlooked a couple of sum operations. The original expression is a triple sum over matrices of dimension 2 or 3, so the result should be a scalar rather than 2 D array.

Sign in to comment.


John BG
John BG on 27 Feb 2017
Edited: John BG on 27 Feb 2017
John
with the additional details perhaps the following would qualify to accepted answer
1.
for 1 product
B=[2 1;3 2;1 4]; % miles between factories and distributors
% B(i,j) miles between i-th factory and j-th distributor
C=[5e3 2e3;6e3 2e3;1e3 4.5e3]; % cost per mile between i-th factory and j-th distributor and per batch
cost_distribution_unit=B.*C % distribution costs per unit between i-th factory and j-th distributor
beta=zeros(size(B));
[szc1 szc2]=size(C);
Q=zeros(size(B));
Q_dcp=[200 300;250 150;100 200];
beta_cell={} % calculating all possible combinations of beta
for k=1:1:numel(C)
L=combinator(6,k,'c');
[szL1 szL2]=size(L);
for s=1:1:szL1
beta_cell=[beta_cell; L(s,:)];
end
end
% beta_cell now contains all possible factory-distributor connection combinations
total_cost=zeros(numel(beta_cell),1);
total_cost_unit=zeros(numel(beta_cell),1);
for k=1:1:numel(beta_cell)
% [ni nj]=ind2sub(size(B),beta_cell{k}) % ni: what factory nj: to what distributor, not used
beta=zeros(size(B));
beta(beta_cell{k})=1 % generate k-th mask
Qtemp=Q_dcp.*beta
cost_distribution_unit_ij=B.*C.*beta
% cost_unit_mile=B.*beta, not used
Amounts_produced=sum(Qtemp,2)
Amounts_in_distributor=sum(Qtemp)
isequal(sum(Amounts_produced),sum(Amounts_in_distributor)) % check Amounts_in_distributor==Amounts_produced, else warn cargo loss, not used
total_amount_ij=sum(sum(Qtemp))
total_cost_distribution_unit_ij=sum(sum(B.*C.*beta)) % cost per unit matrix for combination beta_cell{k}
total_cost_distribution_ij=sum(sum(Qtemp.*total_cost_distribution_unit_ij))
total_cost(k)=total_cost_distribution_ij
cost_unit_ij=total_cost_distribution_ij/total_amount_ij
total_cost_unit(k)=cost_unit_ij
cost_link_ij = Qtemp.*cost_distribution_unit_ij
cost_link_distributor_ij=sum(cost_link_ij)
cost_link_factory_ij=sum(cost_link_ij,2)
end
figure(1);stem(total_cost);title('total cost - 1 product')
.
.
figure(2);stem(total_cost_unit);title('total cost/unit - 1 product')
.
.
2.
3 products
Q_dcp_product_1=[200 300;250 150;100 200];
Q_dcp_product_2=[250 150;100 200;200 300];
Q_dcp_product_3=[100 200;200 300;250 150];
N=3;
Q_dcp_3products=zeros([size(Q_dcp_product_1) N])
Q_dcp_3products(:,:,1)=Q_dcp_product_1;
Q_dcp_3products(:,:,2)=Q_dcp_product_2;
Q_dcp_3products(:,:,3)=Q_dcp_product_3;
% using same beta combinations on all products
total_cost_Nproducts=zeros(N,numel(beta_cell));
total_cost_unit_Nproducts=zeros(N,numel(beta_cell));
for n=1:1:N
for k=1:1:numel(beta_cell)
beta=zeros(size(B));beta(beta_cell{k})=1 % generate k-th beta mask
Q_dcp=Q_dcp_3products(:,:,n);
Qtemp=Q_dcp.*beta
cost_distribution_unit_ij=B.*C.*beta
total_amount_ij=sum(sum(Qtemp))
total_cost_distribution_unit_ij=sum(sum(B.*C.*beta)) % cost per unit matrix for combination beta_cell{k}
total_cost_distribution_ij=sum(sum(Qtemp.*total_cost_distribution_unit_ij))
total_cost_Nproducts(n,k)=total_cost_distribution_ij % update total cost
cost_unit_ij=total_cost_distribution_ij/total_amount_ij
total_cost_unit_Nproducts(n,k)=cost_unit_ij % update cost/unit
% cost_link_ij = Qtemp.*cost_distribution_unit_ij, not used
% cost_link_distributor_ij=sum(cost_link_ij)
% cost_link_factory_ij=sum(cost_link_ij,2)
end
end
3.
combining all products
the total distribution cost, for each beta:
sum(total_cost_Nproducts)
and the total distribution cost/unit, for each beta:
sum(total_cost_unit_Nproducts)
now you could proceed selecting all those elements of beta_cell that only have for instance 4 non zero elements and compare, or find out whether it pays off to have the plants alternating shipments to distributors throughout the week.
if you find this answer useful would you please be so kind to mark my answer as Accepted Answer?
To any other reader, please if you find this answer of any help solving your question,
please click on the thumbs-up vote link,
thanks in advance
John BG

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!