基于异构综合学习粒子群优化波能转换器位置附matlab代码

% Based on the article by Wu (1995) % 'The interaction of water waves with a group of submerged spheres' % All variables have the same notations as in the paper % Infinite water depth % Waves come from right % Calculation of hydrodynamic coefficients and excitation force % 04/06/2015 Algorithms is changed to solve everything in one big matrix (now coupled coefficients are calculated in a right way) % Based on \Debugging\ArraySubmerged_v9.m % 12/06/2015 Changed boundaries for the inetgrals evaluation in quadgk % 19/06/2015 Modified for parallel computation, now frequency and % wavenumber are not in the wave structure, but separate variables % 31/08/2015 Added varargin % - varargin{1} - flag for the form of the output variables 0 - (l, i, j, k), 1 - as (numSphere*3 x numSphere*3) matrix % 11/09/2015 Corrected the calculation of the excitation force % 11/07/2016 Added the range of wave angles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % OUTPUT, Depends on the varargin{1} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Flag = 0 % Coefficients of the sphere l in mode i due to the motion of sphere k in mode j % addedMass = myu(l,i,k,j) % dampingCoef = lambda(l,i,k,j) % excForce = F_exc(l,j,M) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Flag = 1, N - numSphere*3, M - numWaveAngle % addedMass = A(N, N) % dampingCoef = B(N, N) % excForce = F(N, M) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [addedMass, dampingCoef, excForce] = arraySubmergedSphereParfor(array, wave, frequency, waveNumber, numApprox, varargin) if nargin == 5 flag2D = 0; else flag2D = varargin{1}; end numSphere = array.number; ro = wave.waterDensity; % K = wave.waveNumber; % w = wave.frequency; K = waveNumber; w = frequency; alpha = wave.angle + pi; numWaveAngle = length(alpha); radiusSphere = array.radius; xSphere = array.sphereCoordinate(1,:); ySphere = array.sphereCoordinate(2,:); zSphere = array.sphereCoordinate(3,:); %disp(['***position=',mat2str(round(array.sphereCoordinate))]) num_nm = (1+numApprox)/2*numApprox; num_lnm = numSphere*num_nm; % Matrices identification coef_AB = zeros(2*num_lnm); fh = zeros(2*num_lnm, numSphere, 3); fhX = zeros(2*num_lnm, numWaveAngle); myu = zeros(numSphere, 3, numSphere, 3); lambda = zeros(numSphere, 3, numSphere, 3); F_exc = zeros(numSphere, 3, numWaveAngle); M = containers.Map('KeyType','double','ValueType','any'); fctrM = containers.Map('KeyType','double','ValueType','double'); lgM = containers.Map('KeyType','double','ValueType','any'); for l = 1:numSphere sphere_l = [xSphere(l) ySphere(l) zSphere(l)]'; a_l = radiusSphere(l); z_l = zSphere(l); y_l = ySphere(l); x_l = xSphere(l); exp_t = exp(K*z_l+1i*K*(x_l*cos(alpha)+y_l*sin(alpha))); for n = 0:numApprox-1 for m = 0:n % After Eq.(12) if m == 0 eps_m = 1; else eps_m = 2; end % Index (row) of the element for coefficients in front of A and B ind_lnm1 = num_nm*(l-1)+(1+n)/2*n+m+1; % Eq.(19a) ind_lnm2 = num_lnm+num_nm*(l-1)+(1+n)/2*n+m+1; % Eq.(19b) coef_AB(ind_lnm1,ind_lnm1) = coef_AB(ind_lnm1,ind_lnm1) -(n+1)/a_l; coef_AB(ind_lnm2,ind_lnm2) = coef_AB(ind_lnm2,ind_lnm2) -(n+1)/a_l; for k = 1:numSphere % Eq.(20) fh(ind_lnm1,k,1) = -(n==1)*(m==1)*(k==l); fh(ind_lnm1,k,2) = 0; fh(ind_lnm1,k,3) = (n==1)*(m==0)*(k==l); fh(ind_lnm2,k,1) = 0; fh(ind_lnm2,k,2) = -(n==1)*(m==1)*(k==l); fh(ind_lnm2,k,3) = 0; end % Eq.(29a)-(29b) if isKey(fctrM,n+m) fnpm = fctrM(n+m); else fnpm = factorial(n+m); fctrM(n+m) = fnpm; end fhX(ind_lnm1,:) = n*eps_m*w*(-1i)^(m+1)*exp_t.*cos(m*alpha)*(K*a_l)^(n-1)/fnpm; % 11/07/2016 NYS fhX(ind_lnm2,:) = n*eps_m*w*(-1i)^(m+1)*exp_t.*sin(m*alpha)*(K*a_l)^(n-1)/fnpm; % 11/07/2016 NYS for lam = 1:numSphere a_lam = radiusSphere(lam); sphere_lam = [xSphere(lam) ySphere(lam) zSphere(lam)]'; vect = sphere_l - sphere_lam; % Eq.(17) z_lam = zSphere(lam); y_lam = ySphere(lam); x_lam = xSphere(lam); r_laml = norm(vect); Rlaml = sqrt((x_l-x_lam)^2+(y_l-y_lam)^2); theta_laml = acos(vect(3)/r_laml); beta_laml = atan2(vect(2),vect(1)); for nd = 0:numApprox-1 % n dash % cached calculation ind_lamnd = num_nm*(lam-1)+(1+nd)/2*nd; coef_d = eps_m*n*a_lam^(nd+1)*a_l^(n-1); for md = 0:nd % m dash % Index (column) of the element for coefficients in front of A and B ind_lamndmdA = ind_lamnd+md+1; ind_lamndmdB = num_lnm+ind_lamnd+md+1; mdmm=md-m; mdpm=md+m; % hashcode for the key mparas = 173*(173*(173*(23*(n+nd)) + (z_l+z_lam)) + mdmm) + Rlaml; pparas = 173*(173*(173*(23*(n+nd)) + (z_l+z_lam)) + mdpm) + Rlaml; % funM = @(x)x.^(n+nd).*(x+K)./(x-K).*exp(x*(z_l+z_lam)).*besselj(m-md, x*Rlaml); % funP = @(x)x.^(n+nd).*(x+K)./(x-K).*exp(x*(z_l+z_lam)).*besselj(mdpm, x*Rlaml); if isKey(M, mparas) int_mMmd = M(mparas); else funM = @(x)x.^(n+nd).*(x+K)./(x-K).*exp(x*(z_l+z_lam)).*besselj(mdmm, x*Rlaml); int_mMmd1 = quadgk(funM, 0, K*1.5, 'Waypoints', [K*0.7, K+0.1*K*1i, K*1.25]); int_mMmd2 = integral(funM, K*1.5, inf); % changed 12/06/2015 int_mMmd = int_mMmd1 + int_mMmd2; % int_mMmd = apxIntegral(funM, 0, inf, 'Waypoints', [K*0.7, K+0.1*K*1i, K*1.25]); % changed 12/06/2015 M(mparas) = int_mMmd; end if (mdmm) == (mdpm) int_mPmd = int_mMmd; else if isKey(M, pparas) int_mPmd = M(pparas); else funP = @(x)x.^(n+nd).*(x+K)./(x-K).*exp(x*(z_l+z_lam)).*besselj(mdpm, x*Rlaml); int_mPmd1 = quadgk(funP, 0, K*1.5, 'Waypoints', [K*0.7, K+0.1*K*1i, K*1.25]); int_mPmd2 = integral(funP, K*1.5, inf); int_mPmd = int_mPmd1 + int_mPmd2; % int_mPmd = apxIntegral(funP, 0, inf, 'Waypoints', [K*0.7, K+0.1*K*1i, K*1.25]); % changed 12/06/2015 M(pparas)=int_mPmd; end end % I11 = (-1)^m*pi*( cos((m-md)*beta_laml)*int_mMmd+(-1)^md*cos((mdpm)*beta_laml)*int_mPmd); % I12 = (-1)^m*pi*(-sin((m-md)*beta_laml)*int_mMmd+(-1)^md*sin((mdpm)*beta_laml)*int_mPmd); % I21 = (-1)^m*pi*( sin((m-md)*beta_laml)*int_mMmd+(-1)^md*sin((mdpm)*beta_laml)*int_mPmd); % I22 = (-1)^m*pi*( cos((m-md)*b