Implementing a sum with summands of unequal spacing inside an integral

2 views (last 30 days)
Dear all,
for a university study, I need to implement the following equations which gives me some headaches as I'm not a MATLAB expert.
For reasons of good scientific practice, it shall be noted that the following equations are citated from: "Analytical calculation of D- and Q-axis Inductance for interior Permanent Magnet Motors Based on Winding function Theory", Liang et al. in Energies 2016, 9, pp.580-591.
where:
N, k_wv, p, I, Rs, Rr and L, are constants. Same holds for U_d1, U_d2, alpha, beta and mu_0 (see below). Besides, B_d is a combination of multiple Fourier-series:
There are three big issues that I'm struggeling with:
1) While j=1,...,n, for v holds v=1,5,7,11. Therefore I wasn't able to use the "symsum" command as the definition of lower and upper bound for v would include that every integer within 1-11 would be considered.
2) I'm not sure if I can just say (1:j) meaning j sould go from 1 to J and define J=100 later on because...
3) The error message I'm receiving is " Index exceeds number of array elements (1)" in the line "DataStorage(1,i)=B_d...
Although I'm already searching for a while to solve 3), I'm almost sure that my approach isn't purposeful. So any help is highly appreciated!
Finally, please notice that I also would highly appreciate every piece of advice considering only a part of this mess here (e.g. concerning the "v-problem"). Thanks in advance.
This is what I tried:
%% Parameter
p=4;N=32;I=1;alpha=0.742;beta=0.886;Rs=0.091;Rr=0.090;mu_0=4*pi*10^(-7);
L=0.15;k_w=0.95;W=[1,5,7,11];U_d1=6.0368;U_d2=3.3872;
%% Calculation
syms theta j
B_d=@(v,j)(3*N*k_w*I/(v*p*pi)*cos(v*theta)-U_d1*4/pi*sum(sin((1:j)*alpha*pi/2)/(1:j)*cos((1:j)*theta))-U_d2*4/pi*sum((sin((1:j)*pi*beta/2)-sin((1:j)*alpha*pi/2))/(1:j)*cos((1:j)*theta)))*mu_0/(Rs-Rr);
summand=@(v,theta)sum(2*N*k_w/(p*pi)*cos(v.*theta)*B_d*(Rs+Rr)/2*L);
for i=1:length(W) %Solve the sum over v and j in the expression for B_d; that should bring a function of theta only.
v=W(i);
DataStorage(1,i)=B_d(v(i),100);
end
B_d=sum(DataStorage, 'all'); %Does this work even there's still a variable in B_d (theta)?
for i=1:length(W)
v=W(i);
L_md=integral(summand(v,theta),0,2*pi)./I;
end
  2 Comments
David Goodmanson
David Goodmanson on 4 Dec 2019
Hi Richard,
Assuming there are no missing parentheses, the expression for Bd contains three separate summations with integer indices; and the summation for nu contains only the four values 1,5,7,11, is that all correct? Is a numerical solution all right as opposed to a symbolic solution (which may well be no simpler than the code for the numerical one?)
Richard
Richard on 5 Dec 2019
Hi David,
thanks for your response. I double-checked on the parentheses and they are correct. Also your statement referring to nu is correct. A numerical solution is definitely all right. In what follows, I may show my modified MATLAB code to you. I think I finally made my way through somehow. Nevertheless, I think there are much more better approaches than mine and I strive to improve my MATLAB skills so I would be grateful if you (or someone else) might show me how a professional would do!
%% Parameter
p=4;N=32;I=1;alpha=0.742;beta=0.886;Rs=0.091;Rr=0.090;mu_0=4.*pi.*10^(-7);
L=0.15;k_w=0.95;W=[1,5,7,11];U_d1=6.0368;U_d2=3.3872;
J=100;
%% Calculation
syms theta
Fcos=@(v,theta)sum((3.*N.*k_w.*I/(v.*p.*pi)).*cos(v.*theta));
A=@(theta,j)(-U_d1.*4./pi.*sum(sin((1:j).*alpha.*pi/2)./(1:j).*cos((1:j).*theta))-U_d2.*4./pi.*sum((sin((1:j).*pi.*beta./2)-sin((1:j).*alpha.*pi/2))./(1:j).*cos((1:j).*theta)));
for i=1:length(W) %Solve the sum over v and j in the expression for B_d; that should bring a function of theta only.
v=W(i);
DataStorage(i)=([Fcos(v,theta)]);
end
SumFcos=(sum(DataStorage,'all'));
A=(A(theta,J));
B_d=(SumFcos+A).*mu_0./(Rs-Rr);
summand=@(v,theta)(2.*N.*k_w./(p.*pi).*cos(v.*theta));
for i=1:length(W)
v=W(i);
Summand(i)=[summand(v,theta)];
end
Summand=sum(Summand,'all').*B_d.*(Rs+Rr)./2.*L;
L_md=eval(int(Summand,theta,0,2.*pi)./I);

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 5 Dec 2019
Edited: David Goodmanson on 5 Dec 2019
Hi Richard,
I don't claim to be a professional although I did use Matlab at work, and here is how I would do this. The main thing has nothing to do with Matlab. It's just to recognize that the integral is from 0 to 2 pi and contains cos(nu*theta) in the integrand and another cos(j*theta) or cos(nu*theta) [separate sum so this is an independent nu from the first one] in Bd. Therefore orthogonality applies, so all that survives are the factors of cos(nu*theta) in the integrand times the factors of cos(nu*theta) in Bd, both for the same nu. Of the sum j=1 to 100, only j = [1 5 7 11] matter.
Integration 0 to 2pi gives an extra factor of pi since the average value of cos^2 = 1/2. The following result is done numerically and agrees with yours (caveat below).
p=4;N=32;I=1;alpha=0.742;beta=0.886;Rs=0.091;Rr=0.090;mu_0=4.*pi.*10^(-7);
L=0.15;k_w=0.95;W=[1,5,7,11];U_d1=6.0368;U_d2=3.3872;
nu = W; % notation
% include everything except cos(nu*theta) or cos(j*theta) factors
intgr = (2*N*(k_w/(p.*pi))*((Rr+Rs)/2)*L/I)./nu;
% toss in overall constant for Bd later on
Bd = 3*N*(k_w/(p.*pi))*I./nu -(4*U_d1/pi)*sin(nu*alpha*pi/2)./nu ...
-(4*U_d2/pi)*(sin(nu*beta*pi/2)-sin(nu*alpha*pi/2))./nu;
result = pi*sum(intgr.*Bd*(mu_0/(Rs-Rr)))
The caveat is, you seem to be missing a factor of 1/nu in the integrand expression. Insert that, and the methods agree.
  1 Comment
Richard
Richard on 5 Dec 2019
Marvelous! I admit that I first had to look up about the orthogonality of the cos(.)-expressions but of course your're totally right. I think it really shows that one (me :) )should rather look at problems like that with "eyes open" than just start to code the equations!
Anyway, sincere thanks for your efford and for highlighting the missing nu to me!
Best regards
Richard

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!