# Optimization of a Variable in a Loop

5 views (last 30 days)
Peter Vella on 1 Sep 2019
Answered: Thomas Barbier on 2 Oct 2019
I have the below script as part of a larger script to calculate expander power over a yearly data set of available stroage of comrpessed air.
I need to change "m" from a constant (20.2164) such that I obtain aspecific value of "Pwr_Bk1_St2" (summed up across r).
What would be the ebst way for this? Do I use linear/ non linear optimisation? Or do I need to include if-then loops?
Thanks
+++
ts=7; % Time step for calculations (>1).
Pi=P1e(7,2); PRi=PR(PRr,2); a1=(Pi-(Pi/PRi))/ts; eff_c=.75; k=1.4; b=1; Cp=1.005; m=20.2164; % Setting up the "parameters", i.e. the "constants".
P=zeros(ts,2); % Initialising the 2D array, P.
% Loop to set up the "pressures" at the established 10-min time steps.
n=((k/(k-1))*eff_c)/((k/(k-1))*eff_c-1);
for r=1:1:ts
b=b-1;
for c=1:2
P2e(r,c)=(Pi-b*a1);
b=b+1;
T12(1)=(330+273);
T2a(r)=((P2e(r,2)/P2e(r,1))^((1.27-1)/1.27))*T12(r);
if r<=6
T12(r+1)=T2a(r);
end
T12(1)=(330+273);
T2s2(r)=((P2e(r,2)/P2e(r,1))^((1.4-1)/1.4))*T12(r);
Pwr_Bk1_St2(r)=(m/7)*Cp*(T12(r)-T2s2(r));
end
end

Walter Roberson on 1 Sep 2019
Generally speaking,
desired_m = fzero(@(m) some_function(m) - Pwr_Bk1_St2_you_want, initial_value_of_m )
Peter Vella on 1 Sep 2019
Thanks Walter.
Have been learning MATLAB for only 3 months .. Will give it a try.

Thomas Barbier on 2 Oct 2019
Hi,
I would personnaly use the lsqnonlin function from the optimisation toolbox, as your problem seem quite non-linear.
m0 = 20.2164;
m = lsqnonlin(myFunction - desired_Pwr_Bk1_St2, m0);
function Pwr_Bk1_St2 = myFunction(m)
ts=7; % Time step for calculations (>1).
Pi=P1e(7,2); PRi=PR(PRr,2); a1=(Pi-(Pi/PRi))/ts; eff_c=.75; k=1.4; b=1; Cp=1.005; m=20.2164; % Setting up the "parameters", i.e. the "constants".
P=zeros(ts,2); % Initialising the 2D array, P.
% Loop to set up the "pressures" at the established 10-min time steps.
n=((k/(k-1))*eff_c)/((k/(k-1))*eff_c-1);
for r=1:1:ts
b=b-1;
for c=1:2
P2e(r,c)=(Pi-b*a1);
b=b+1;
T12(1)=(330+273);
T2a(r)=((P2e(r,2)/P2e(r,1))^((1.27-1)/1.27))*T12(r);
if r<=6
T12(r+1)=T2a(r);
end
T12(1)=(330+273);
T2s2(r)=((P2e(r,2)/P2e(r,1))^((1.4-1)/1.4))*T12(r);
Pwr_Bk1_St2(r)=(m/7)*Cp*(T12(r)-T2s2(r));
end
end
Pwr_Bk1_St2 = sum(Pwr_Bk1_St2);
end

Translated by