How do I write the code for open luup nd closed luup ?

To simulate the open loop and closed loop dynamic behavior of STHE

Answer 1

[Pablo Escobar]

The answer for open luup is :-

clc
close all
clear all
h = 0.5; %in minutes
run_count = 201; %time count

%FIXED PARAMETERS
cp = 1; %cal/gm C
rho = 1e6; %gm/cubic meters
cpc = 1; %cal/gm C
rhoc = 1e6; %gm/cubic meters
a = 1.41e5;
b = 0.5;
Tcs = 25;

%VARYING PARAMETERS
v = 3.5; %cubic meters
fs = 0.085; %cubic meters/min
fcs = 0.44; %cubic meters/min
Ths= 150;

alp = fs/v;
beta1 = a*fcs^(b+1)/(rho*cp*v);
beta2 = a*fcs^b/(2*rhoc*cpc);
beta = beta1/(fcs + beta2);
Ts = (alp*Ths + beta*Tcs)/(alp + beta);

Tc = Tcs; %constant
Th = Ths*ones(run_count,1); %disturbance
f = fs*ones(run_count,1); %disturbance
fc = fcs*ones(run_count,1); %manipulated variable
t = zeros(run_count,1);

T = Ts*ones(run_count,1);
noise_std = 0.07;
Tm = zeros(run_count,1);
Tm(1) = T(1) + noise_std*randn;

%introduction of disturbance
n = input('Specify the variable for introducing step change. \n Enter 1 for coolant flow rate. \n Enter 2 for inlet stream temperature.\n');
step_time = 21;
if n == 1
    delta_fc = 0.1*fcs; % introducing +10% change in flow rate
    fc(step_time:end) = delta_fc + fc(step_time:end);
else
    delta_t = +10; %introducing a +10 deg. C change in inlet temperature
    Th(step_time:end) = Th(step_time:end) + delta_t;
end

for k = 1:run_count-1
    alp_ = f(k)/v;
    beta1_ = a*fc(k)^(b+1)/(rho*cp*v);
    beta2_ = a*fc(k)^b/(2*rhoc*cpc);
    beta_ = beta1_/(fc(k) + beta2_);    
    T(k+1) = T(k) + h*[alp_*(Th(k) - T(k)) - beta_*(T(k) - Tc)];    
    Tm(k+1) = T(k+1) + noise_std*randn;
    t(k+1)=k*0.5;
end

set(0,'DefaultLineLineWidth',2)
set(0,'DefaultaxesLineWidth',2)
set(0,'DefaultaxesFontSize',16)
set(0,'DefaultTextFontSize',16)
set(0,'DefaultaxesFontName','arial')

figure(1), plot(t, T,'b', t, Tm, 'g',t, Ts*ones(run_count,1),'r'), grid;
xlabel('Time (min)')
ylabel('STHE Temperature (deg C)')
title('Controlled Output (T)')
legend('True Temperature','Measured Temperature','Initial S.S. Temperature')

figure(2), stairs(t, fc), grid;
xlabel('Time (min)')
ylabel('Coolant Flow (cubic meters/min)')
title('Manipulated Input (F_c)')

figure(3), stairs(t, Th), grid;
xlabel('Time (min)')
ylabel('Inlet Stream Temperature (deg C)')
title('Disturbance (T_h)')

diff = [mean(Tm(end-10:end))-mean(Tm(1:10))];
delta = 0.63*diff*ones(run_count,1);
figure(4), plot(t, Tm - Ts,'b' ,t, delta,'r'), grid;
xlabel('Time (min)')
ylabel('Deviation STHE Temperature (deg C)')
title('Deviation Controlled Output (T)')
legend('Measured Output','63.2% of Change in Output')

figure(5), stairs(t, fc-fcs), grid;
xlabel('Time (min)')
ylabel('Deviation Coolant Flow (cubic meters/min)')
title('Deviation Manipulated Input (Fc)')

figure(6), stairs(t, Th - Ths), grid;
xlabel('Time (min)')
ylabel('Deviation Inlet Temperature (deg C)')
title('Deviation Disturbance (T)')

clc

if n == 1
    K_mani = diff/delta_fc

else
    K_dist = diff/delta_t

end

Answer 2

[Pablo Escobar]

The answer for closed luup is :-

clear all
clc
close all

% Initia;lizing of variables

% set sampling time / integration interval & simulation runtime
int_T = 0.5; % Integration interval/sampling time = 0.5 min
run_count = 201; % Simulation time count

%Initializing of Model parameters

STHE.V = 2.1; % Volume of STHE (m3)
STHE.Cp = 1.0; % Cp vale (cal/g-C)
STHE.rho = 10^06; % Density (g/m3)
STHE.Cpc = 1.0; % Cp vale of coolant (cal/g-C)
STHE.rhoc = 10^06; % Density of coolant (g/m3)
STHE.a = 1.41*(10^5); % Heat transfer eqn parameter (cal/min-C)
STHE.b = 0.5; % Heat transfer eqn parameter

%Define steady state input values at desired operating point

STHE.Ths = 150; % Hot fluid inlet temperature (C)
STHE.Fs = 0.085; % Hot fluid inlet flow rate
STHE.Tcins = 25; % Coolant temperature (C)
STHE.Fcs = 0.5; % Coolant flow rate 

%Find steady state value of state varaible/ controled  output

alfa = STHE.Fs/STHE.V ;
num = STHE.a*STHE.Fcs^(STHE.b+1)/(STHE.rho*STHE.Cp*STHE.V);
Den = STHE.a*STHE.Fcs^STHE.b/(2*STHE.rhoc*STHE.Cpc);
beta = num /(STHE.Fcs + Den);
STHE.Ts = (alfa*STHE.Ths + beta*STHE.Tcins)/(alfa+beta);

% define initial state for dynamic simulation

STHE.Tcin = STHE.Tcins; %initial coolant temp(assumed constant)
STHE.F = STHE.Fs*ones(run_count,1); %Inlet hot fluid flow rate(disturbabnce)
STHE.Th = STHE.Ths*ones(run_count,1); %hot fluid temp(dist)
STHE.Fc = STHE.Fcs*ones(run_count,1); %initial cooant flow(manipulated i/p)

Time(1) = 0;
STHE.T(1) = STHE.Ts; %Initliaize the 1st element in temp array

noise_std=0.07; %Standard deviation of measurement noise
STHE.Tm(1) = STHE.T(1) + noise_std*randn; %1st temp measurement

% PID Controller Parameters

PID.Kc = -1/20 ; % From direct synthesis method (from open loop expt)
% PID.Kc = input('Enter the Kc value obtained from open loop simulation : \n');
PID.Ti = 10 ; % From direct synthesis method  (from open loop expt)
% PID.Ti = input('Enter the Ti value obtained from open loop simulation : \n');
PID.error = zeros(run_count,1);
PID.ip_max = 2 * STHE.Fcs; % Upper bound on manipulated input
PID.ip_min = 0.0001 ; % Sure to put 0??

fprintf('\n\n\t Closed loop Dynamic Simulation \n\n\t')
change_type = input('Specify your choice :: \n PRESS 1 for Servo case \n PRESS 2 for Regulatory case \n\n' ) ;
setpoint = STHE.Ts * ones(run_count,1) ;

if ( change_type ==1 )
delta_setpt = 5 ; % Introduce setp change (SERVO)
step_time = 21 ; 
setpoint(step_time:end) = setpoint(step_time:end) - delta_setpt;
elseif (change_type ==2)
delta_Th = 10; % Introduce step change (REGULATORY)
step_time = 21 ;
STHE.Th(step_time:end) = STHE.Th(step_time:end) + delta_Th ;
else
    fprintf('This is not a correct choice, please choose between 1 and 2 \n');
    break 
end

% Closed Dynamic Simulation Loop

for k =1 :run_count-1

% Controller calculations at instant t = k*T
error(k) = setpoint(k) - STHE.Tm(k) ; % Calculate controlled error

if ( k>1 )
% PI Controller Calculations
P_contribution = PID.Kc * (error(k) - error(k-1) ) ;
% I_contribution = (PID.Kc*int_T / PID.Ti) * error(k) ;
I_contribution = 0;
delta_Fc = P_contribution + I_contribution ;
STHE.Fc(k) = STHE.Fc(k-1) + delta_Fc ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The IF-ELSE loop is provided to ensure that the measured value stick in
% the limited span (between the MAX and the MIN values)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (STHE.Fc(k) > PID.ip_max)
STHE.Fc(k) = PID.ip_max ;
elseif (STHE.Fc(k) < PID.ip_min)
STHE.Fc(k) = PID.ip_min ;
end
end

      %Process simulation from time t=k*t to t =(k+1)*T using explicit eulter
      %integration

      Time(k+1)=k*int_T;

      alfa=STHE.F(k)/STHE.V;
      num = STHE.a*STHE.Fc(k)^(STHE.b+1)/(STHE.rho*STHE.Cp*STHE.V);
      Den = STHE.a*STHE.Fc(k)^STHE.b/(2*STHE.rhoc*STHE.Cpc);
      beta = num/(STHE.Fc(k)+Den);
      %Compute Derivation at t=time(k)
      STHE.Q=beta*(STHE.T(k) - STHE.Tcin);
      dT_by_dt= alfa*(STHE.Th(k) - STHE.T(k))-STHE.Q;
      STHE.T(k+1) = STHE.T(k) + int_T*dT_by_dt; %Euler integration

      %Generate teperature measurement corrupted with noise
      STHE.Tm(k+1) = STHE.T(k+1)+noise_std*randn;
  end
  % Display simulation results

%Initalize plotting style
set(0,'DefaultLineLineWidth',2)
set(0,'DefaultaxesLineWidth',2)
set(0,'DefaultaxesFontSize',16)
set(0,'DefaultTextFontSize',16)
set(0,'DefaultaxesFontName','arial')

% DISPLAY SIMULATION RESULTS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% DISPLAYING RESULTS IN ABSOLUTE VALUES

figure(1)
plot(Time,setpoint,'--',Time,STHE.T,'-',Time,STHE.Tm,'-.'),grid
xlabel('Time(min)')
ylabel('STHE temperature (deg C)')
title('Controlled output (T)')
legend('Setpoint','True output','Measured output')

figure(2)
stairs(Time,STHE.Fc),grid
xlabel('Time(min)')
ylabel('Coolant flow (m3/s)')
title('Manipulated input (F_c)')

figure(3)
stairs(Time,STHE.Th),grid
xlabel('Time(min)')
ylabel('Hot fluid inlet temperature (deg C)')
title('Disturbance (T_h)')

% DISPLAYING RESULTS IN PERTURBATION VALUES

figure(4)
plot(Time,setpoint - STHE.Ts,'--',Time,STHE.T - STHE.Ts,'-',Time,STHE.Tm - STHE.Ts,'-.'),grid
xlabel('Time(min)')
ylabel('Deviation in STHE temperature (deg C)')
title('Deviation in Controlled output (T)')
legend('Setpoint','True output','Measured output')

figure(5)
stairs(Time,STHE.Fc - STHE.Fcs),grid
xlabel('Time(min)')
ylabel('Deviation in Coolant flow (m3/s)')
title('Deviation in Manipulated input')

figure(6)
stairs(Time,STHE.Th - STHE.Ths),grid
xlabel('Time(min)')
ylabel('Deviation in Hot fluid inlet temperature (deg C)')
title('Deviation in Disturbance')

%%%%%%%
% END %
%%%%%%%

Answer 3

[Pablo Escobar]

PFA Answers