Please assist to resolve "magaa and mag2a" codes problems; probably FUN and fsolve, not properly sets. The mag2 has to generate 12 (Qy-variable) questions by for-end. Thanks

조회 수: 4 (최근 30일)
I have two file execution files, MAGa is main and mag2a is responsible for produce and generate 12 non-linear equations, and probably could goes to 100. The MAGa file has to call mag2a and produce the output in table format. I have failed to known what is main problem, include set up FUN in mag2 file and solver (non-linear) I has proposed. See the attached files and executed output; appreciate your assistance
FILE MAGa and mag2a
MAGa
function [Qy, FUN] = MAGa(Qy0)
tic, n=12; pipes=n; dia=n; LP=4; cCc=7; HLD=-0; ee=0.04; vs=1.00001*10^-6; format shortEng
D=[1.120, 0.059, 0.019, 0.033, 0.060, 0.052, 0.023, 0.029, 0.017, 0.018, 0.058, 0.150]; L=[50, 50, 40, 50, 40, 50, 60, 60, 50, 50, 60, 40];
FUN = @mag2a; %
options = optimoptions('fsolve','Display','iter',MaxIterations=10000000, MaxFunctionEvaluations=10000000, Algorithm='levenberg-marquardt');
% Qy0(1:12,1)=0.30;
[Qy]=fsolve(FUN, Qy0, options);
Len=L(1,:)'; Dia=D(1,:)'; vel=4*Qy(:,1)./(3.14.*Dia(:,1).^2); Re=4.*Qy(:,1)./(3.14*vs.*Dia(:,1));
fy11=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)*0.0001.^0.5)); fy1=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)).*fy11.^0.5); fy=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)).*fy1.^0.5); f=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)).*fy.^0.5); ff=(1./f).^2;
Data1=zeros(n,cCc); Data1(:,1)=Len(:,1); Data1(:,2)=Dia(:,1); Data1(:,3)=Qy(:,1); Data1(:,4)=vel(:,1); Data1(:,5)=ff(:,1); Data1(:,6)=Re(:,1); % Data1,
StrNum={'string','double','double','double','double','double','double','double'}; subTitl={'Pipes [x-z]', 'Length [m]', 'Diameter [m]', 'Flow rate [m3/s]', 'Velocity [m/s]', 'Friction, f', 'Reynold No.', 'Headlosses [m]'}; %
tz=table('size', [n, cCc+1], 'variabletypes',StrNum, 'variablenames', subTitl); Pipexz=["1P(N1-N2)","2P(N1-N2)","3P(N1-N2)","4P(N1-N2)","5P(N1-N2)","6P(N1-N2)","7P(N1-N2)","8P(N1-N2)","9P(N1-N2)","10P(N1-N2)","11P(N1-N2)","12P(N1-N2)"];
tz{:,1}=Pipexz'; tz{:,2:cCc+1}=Data1(:,1:cCc), toc % MIDsrt=sortrows(tz,sn); MIDfn=MIDsrt; MIDfn(:,end)=[ ];
end
mag2a
function [FUN] = mag2a(D,L,LP,n,ee,vs,matr,HLD)
matr=[1,0,0,0,1,-1,0,0,0,0,0,-1; 0,1,1,-1,-1,0,0,0,0,0,0,0; 0,0,0,1,0,0,1,-1,-1,0,0,0; 0,0,0,0,0,1,0,1,0,-1,-1,0]; LP=4; n=12;
fprintf('FUN=@(Qy)[');
for j=1:LP, j;
if (j<=4-0)
for k=1:n, k;
if sum(matr(j,1:k))~=sum(matr(j,:))&&matr(j,k)~=0 fprintf('%d*(Qy(%d)*abs(Qy(%d)*8*L(1,%d)*1./(9.81*3.14.^2*D(1,%d).^5)*(-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...\n',matr(j,k),k,k,k,k,k,k,k,k,k,k,k,k,k), end
if sum(matr(j,1:k))==sum(matr(j,:))&&matr(j,k)~=0 fprintf('%d*(Qy(%d)*abs(Qy(%d)*8*L(1,%d)*1./(9.81*3.14.^2*D(1,%d).^5)*(-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2;\n',matr(j,k),k,k,k,k,k,k,k,k,k,k,k,k,k),
end, end, end
end
fprintf('Qy(1)-Qy(2)-Qy(5)-0.03; Qy(2)-Qy(3)-0.04; Qy(3)+Qy(4)-Qy(7)-0.03; Qy(6)+Qy(5)-Qy(4)-Qy(8)-0.05; Qy(7)+Qy(9)-0.04; Qy(8)+Qy(10)-Qy(9)-0.03; Qy(11)-Qy(10)-0.04; Qy(12)-Qy(6)-Qy(11)-0.04];\n'); % 1/(Qy(7)).^0.5+4*log(ee./(3.7*D(1,1))+2.54*3.14*D(1,1)*vs./(4*Qy(1)*(Qy(7)).^0.5));...
end
OUTPUT OF EXECUTION (OUTPUTS)
Qy0(1:12,1)=0.30; MAGa(Qy0)
FUN=@(Qy)[1*(Qy(1)*abs(Qy(1)*8*L(1,1)*1./(9.81*3.14.^2*D(1,1).^5)*(-2*log10(ee./(3.7.*D(1,1))+2.51*3.14*vs.*D(1,1)./(4.*Qy(1,1)*((-2*log10(ee./(3.7.*D(1,1))+2.51*3.14*vs.*D(1,1)./(4.*Qy(1,1)*((-2*log10(ee./(3.7.*D(1,1))+2.51*3.14*vs.*D(1,1)./(4.*Qy(1,1)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
1*(Qy(5)*abs(Qy(5)*8*L(1,5)*1./(9.81*3.14.^2*D(1,5).^5)*(-2*log10(ee./(3.7.*D(1,5))+2.51*3.14*vs.*D(1,5)./(4.*Qy(1,5)*((-2*log10(ee./(3.7.*D(1,5))+2.51*3.14*vs.*D(1,5)./(4.*Qy(1,5)*((-2*log10(ee./(3.7.*D(1,5))+2.51*3.14*vs.*D(1,5)./(4.*Qy(1,5)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(6)*abs(Qy(6)*8*L(1,6)*1./(9.81*3.14.^2*D(1,6).^5)*(-2*log10(ee./(3.7.*D(1,6))+2.51*3.14*vs.*D(1,6)./(4.*Qy(1,6)*((-2*log10(ee./(3.7.*D(1,6))+2.51*3.14*vs.*D(1,6)./(4.*Qy(1,6)*((-2*log10(ee./(3.7.*D(1,6))+2.51*3.14*vs.*D(1,6)./(4.*Qy(1,6)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(12)*abs(Qy(12)*8*L(1,12)*1./(9.81*3.14.^2*D(1,12).^5)*(-2*log10(ee./(3.7.*D(1,12))+2.51*3.14*vs.*D(1,12)./(4.*Qy(1,12)*((-2*log10(ee./(3.7.*D(1,12))+2.51*3.14*vs.*D(1,12)./(4.*Qy(1,12)*((-2*log10(ee./(3.7.*D(1,12))+2.51*3.14*vs.*D(1,12)./(4.*Qy(1,12)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2;
1*(Qy(2)*abs(Qy(2)*8*L(1,2)*1./(9.81*3.14.^2*D(1,2).^5)*(-2*log10(ee./(3.7.*D(1,2))+2.51*3.14*vs.*D(1,2)./(4.*Qy(1,2)*((-2*log10(ee./(3.7.*D(1,2))+2.51*3.14*vs.*D(1,2)./(4.*Qy(1,2)*((-2*log10(ee./(3.7.*D(1,2))+2.51*3.14*vs.*D(1,2)./(4.*Qy(1,2)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
1*(Qy(3)*abs(Qy(3)*8*L(1,3)*1./(9.81*3.14.^2*D(1,3).^5)*(-2*log10(ee./(3.7.*D(1,3))+2.51*3.14*vs.*D(1,3)./(4.*Qy(1,3)*((-2*log10(ee./(3.7.*D(1,3))+2.51*3.14*vs.*D(1,3)./(4.*Qy(1,3)*((-2*log10(ee./(3.7.*D(1,3))+2.51*3.14*vs.*D(1,3)./(4.*Qy(1,3)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(4)*abs(Qy(4)*8*L(1,4)*1./(9.81*3.14.^2*D(1,4).^5)*(-2*log10(ee./(3.7.*D(1,4))+2.51*3.14*vs.*D(1,4)./(4.*Qy(1,4)*((-2*log10(ee./(3.7.*D(1,4))+2.51*3.14*vs.*D(1,4)./(4.*Qy(1,4)*((-2*log10(ee./(3.7.*D(1,4))+2.51*3.14*vs.*D(1,4)./(4.*Qy(1,4)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(5)*abs(Qy(5)*8*L(1,5)*1./(9.81*3.14.^2*D(1,5).^5)*(-2*log10(ee./(3.7.*D(1,5))+2.51*3.14*vs.*D(1,5)./(4.*Qy(1,5)*((-2*log10(ee./(3.7.*D(1,5))+2.51*3.14*vs.*D(1,5)./(4.*Qy(1,5)*((-2*log10(ee./(3.7.*D(1,5))+2.51*3.14*vs.*D(1,5)./(4.*Qy(1,5)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2;
1*(Qy(4)*abs(Qy(4)*8*L(1,4)*1./(9.81*3.14.^2*D(1,4).^5)*(-2*log10(ee./(3.7.*D(1,4))+2.51*3.14*vs.*D(1,4)./(4.*Qy(1,4)*((-2*log10(ee./(3.7.*D(1,4))+2.51*3.14*vs.*D(1,4)./(4.*Qy(1,4)*((-2*log10(ee./(3.7.*D(1,4))+2.51*3.14*vs.*D(1,4)./(4.*Qy(1,4)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
1*(Qy(7)*abs(Qy(7)*8*L(1,7)*1./(9.81*3.14.^2*D(1,7).^5)*(-2*log10(ee./(3.7.*D(1,7))+2.51*3.14*vs.*D(1,7)./(4.*Qy(1,7)*((-2*log10(ee./(3.7.*D(1,7))+2.51*3.14*vs.*D(1,7)./(4.*Qy(1,7)*((-2*log10(ee./(3.7.*D(1,7))+2.51*3.14*vs.*D(1,7)./(4.*Qy(1,7)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(8)*abs(Qy(8)*8*L(1,8)*1./(9.81*3.14.^2*D(1,8).^5)*(-2*log10(ee./(3.7.*D(1,8))+2.51*3.14*vs.*D(1,8)./(4.*Qy(1,8)*((-2*log10(ee./(3.7.*D(1,8))+2.51*3.14*vs.*D(1,8)./(4.*Qy(1,8)*((-2*log10(ee./(3.7.*D(1,8))+2.51*3.14*vs.*D(1,8)./(4.*Qy(1,8)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(9)*abs(Qy(9)*8*L(1,9)*1./(9.81*3.14.^2*D(1,9).^5)*(-2*log10(ee./(3.7.*D(1,9))+2.51*3.14*vs.*D(1,9)./(4.*Qy(1,9)*((-2*log10(ee./(3.7.*D(1,9))+2.51*3.14*vs.*D(1,9)./(4.*Qy(1,9)*((-2*log10(ee./(3.7.*D(1,9))+2.51*3.14*vs.*D(1,9)./(4.*Qy(1,9)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2;
1*(Qy(6)*abs(Qy(6)*8*L(1,6)*1./(9.81*3.14.^2*D(1,6).^5)*(-2*log10(ee./(3.7.*D(1,6))+2.51*3.14*vs.*D(1,6)./(4.*Qy(1,6)*((-2*log10(ee./(3.7.*D(1,6))+2.51*3.14*vs.*D(1,6)./(4.*Qy(1,6)*((-2*log10(ee./(3.7.*D(1,6))+2.51*3.14*vs.*D(1,6)./(4.*Qy(1,6)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
1*(Qy(8)*abs(Qy(8)*8*L(1,8)*1./(9.81*3.14.^2*D(1,8).^5)*(-2*log10(ee./(3.7.*D(1,8))+2.51*3.14*vs.*D(1,8)./(4.*Qy(1,8)*((-2*log10(ee./(3.7.*D(1,8))+2.51*3.14*vs.*D(1,8)./(4.*Qy(1,8)*((-2*log10(ee./(3.7.*D(1,8))+2.51*3.14*vs.*D(1,8)./(4.*Qy(1,8)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(10)*abs(Qy(10)*8*L(1,10)*1./(9.81*3.14.^2*D(1,10).^5)*(-2*log10(ee./(3.7.*D(1,10))+2.51*3.14*vs.*D(1,10)./(4.*Qy(1,10)*((-2*log10(ee./(3.7.*D(1,10))+2.51*3.14*vs.*D(1,10)./(4.*Qy(1,10)*((-2*log10(ee./(3.7.*D(1,10))+2.51*3.14*vs.*D(1,10)./(4.*Qy(1,10)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+...
-1*(Qy(11)*abs(Qy(11)*8*L(1,11)*1./(9.81*3.14.^2*D(1,11).^5)*(-2*log10(ee./(3.7.*D(1,11))+2.51*3.14*vs.*D(1,11)./(4.*Qy(1,11)*((-2*log10(ee./(3.7.*D(1,11))+2.51*3.14*vs.*D(1,11)./(4.*Qy(1,11)*((-2*log10(ee./(3.7.*D(1,11))+2.51*3.14*vs.*D(1,11)./(4.*Qy(1,11)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2;
Qy(1)-Qy(2)-Qy(5)-0.03; Qy(2)-Qy(3)-0.04; Qy(3)+Qy(4)-Qy(7)-0.03; Qy(6)+Qy(5)-Qy(4)-Qy(8)-0.05; Qy(7)+Qy(9)-0.04; Qy(8)+Qy(10)-Qy(9)-0.03; Qy(11)-Qy(10)-0.04; Qy(12)-Qy(6)-Qy(11)-0.04];
Output argument "FUN" (and possibly others) not assigned a value in the execution with "mag2a" function.
Error in fsolve (line 270)
fuser = feval(funfcn{3},x,varargin{:});
Error in MAGa (line 8)
[Qy]=fsolve(FUN, Qy0, options);
Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.
>>

채택된 답변

Walter Roberson
Walter Roberson 2024년 4월 26일
편집: Walter Roberson 2024년 4월 26일
function [Qy, FUN] = MAGa(Qy0)
tic, n=12; pipes=n; dia=n; LP=4; cCc=7; HLD=-0; ee=0.04; vs=1.00001*10^-6; format shortEng
D=[1.120, 0.059, 0.019, 0.033, 0.060, 0.052, 0.023, 0.029, 0.017, 0.018, 0.058, 0.150]; L=[50, 50, 40, 50, 40, 50, 60, 60, 50, 50, 60, 40];
FUN = '';
eval(mag2a(D,L,LP,n,ee,vs,matr,HLD)); %assigns to FUN
options = optimoptions('fsolve','Display','iter',MaxIterations=10000000, MaxFunctionEvaluations=10000000, Algorithm='levenberg-marquardt');
% Qy0(1:12,1)=0.30;
[Qy]=fsolve(FUN, Qy0, options);
Len=L(1,:)'; Dia=D(1,:)'; vel=4*Qy(:,1)./(3.14.*Dia(:,1).^2); Re=4.*Qy(:,1)./(3.14*vs.*Dia(:,1));
fy11=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)*0.0001.^0.5)); fy1=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)).*fy11.^0.5); fy=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)).*fy1.^0.5); f=-2*log10(ee./(3.7.*Dia(:,1))+2.51*3.14*vs.*Dia(:,1)./(4.*Qy(:,1)).*fy.^0.5); ff=(1./f).^2;
Data1=zeros(n,cCc); Data1(:,1)=Len(:,1); Data1(:,2)=Dia(:,1); Data1(:,3)=Qy(:,1); Data1(:,4)=vel(:,1); Data1(:,5)=ff(:,1); Data1(:,6)=Re(:,1); % Data1,
StrNum={'string','double','double','double','double','double','double','double'}; subTitl={'Pipes [x-z]', 'Length [m]', 'Diameter [m]', 'Flow rate [m3/s]', 'Velocity [m/s]', 'Friction, f', 'Reynold No.', 'Headlosses [m]'}; %
tz=table('size', [n, cCc+1], 'variabletypes',StrNum, 'variablenames', subTitl); Pipexz=["1P(N1-N2)","2P(N1-N2)","3P(N1-N2)","4P(N1-N2)","5P(N1-N2)","6P(N1-N2)","7P(N1-N2)","8P(N1-N2)","9P(N1-N2)","10P(N1-N2)","11P(N1-N2)","12P(N1-N2)"];
tz{:,1}=Pipexz'; tz{:,2:cCc+1}=Data1(:,1:cCc), toc % MIDsrt=sortrows(tz,sn); MIDfn=MIDsrt; MIDfn(:,end)=[ ];
end
function [FUN] = mag2a(D,L,LP,n,ee,vs,matr,HLD)
matr=[1,0,0,0,1,-1,0,0,0,0,0,-1; 0,1,1,-1,-1,0,0,0,0,0,0,0; 0,0,0,1,0,0,1,-1,-1,0,0,0; 0,0,0,0,0,1,0,1,0,-1,-1,0]; LP=4; n=12;
FUN = 'FUN=@(Qy)[';
for j=1:LP;
if (j<=4-0)
for k=1:n, k;
if sum(matr(j,1:k))~=sum(matr(j,:))&&matr(j,k)~=0 FUN = [FUN, sprintf('%d*(Qy(%d)*abs(Qy(%d)*8*L(1,%d)*1./(9.81*3.14.^2*D(1,%d).^5)*(-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+',matr(j,k),k,k,k,k,k,k,k,k,k,k,k,k,k)]; end
if sum(matr(j,1:k))==sum(matr(j,:))&&matr(j,k)~=0 FUN = [FUN, sprintf('%d*(Qy(%d)*abs(Qy(%d)*8*L(1,%d)*1./(9.81*3.14.^2*D(1,%d).^5)*(-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2',matr(j,k),k,k,k,k,k,k,k,k,k,k,k,k,k)]; end
end
end
end
FUN = [FUN, 'Qy(1)-Qy(2)-Qy(5)-0.03; Qy(2)-Qy(3)-0.04; Qy(3)+Qy(4)-Qy(7)-0.03; Qy(6)+Qy(5)-Qy(4)-Qy(8)-0.05; Qy(7)+Qy(9)-0.04; Qy(8)+Qy(10)-Qy(9)-0.03; Qy(11)-Qy(10)-0.04; Qy(12)-Qy(6)-Qy(11)-0.04];']; % 1/(Qy(7)).^0.5+4*log(ee./(3.7*D(1,1))+2.54*3.14*D(1,1)*vs./(4*Qy(1)*(Qy(7)).^0.5));...
end
  댓글 수: 3
Walter Roberson
Walter Roberson 2024년 4월 26일
if sum(matr(j,1:k))~=sum(matr(j,:))&&matr(j,k)~=0 FUN = [FUN, sprintf('%d*(Qy(%d)*abs(Qy(%d)*8*L(1,%d)*1./(9.81*3.14.^2*D(1,%d).^5)*(-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*((-2*log10(ee./(3.7.*D(1,%d))+2.51*3.14*vs.*D(1,%d)./(4.*Qy(1,%d)*(0.0001).^0.5)))).^0.5)))).^0.5))).^-2+',matr(j,k),k,k,k,k,k,k,k,k,k,k,k,k,k)]; end
% 1 2 1 2 3 2 3 2 3 4 3 2 3 4 5 6 54 5 4 5 6 5 67 8 9 A 98 9 8 9 A 9 AB C D E DC D C D E D E D CBA9 8765 432
Your format has two more ( than )
Eng. FM
Eng. FM 2024년 4월 27일
Woo,
You have provided to me unforgotable remark to my analytical task; the complex analysis, I have got new chapter. Be blessed too, Mr Walter.
Probably in few months to come, I will extend the model to handle hundred to thousand nodes and pipes with multiple sources etc
Once again, thanks

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Just for fun에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by