Optimization of a Variable in a Loop

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]

Generally speaking,

desired_m = fzero(@(m) some_function(m) - Pwr_Bk1_St2_you_want, initial_value_of_m )

[Peter Vella]

Thanks Walter.

Have been learning MATLAB for only 3 months .. Will give it a try.

Answer 1

[Thomas Barbier]

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