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

1 次查看（过去 30 天）
Pablo Escobar
Pablo Escobar 2016-9-19
回答： Pablo Escobar 2016-9-19
To simulate the open loop and closed loop dynamic behavior of STHE

请先登录，再回答此问题。

采纳的回答

Pablo Escobar
Pablo Escobar 2016-9-19
编辑：per isakson 2016-9-19
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

更多回答（2 个）

Pablo Escobar
Pablo Escobar 2016-9-19
编辑：per isakson 2016-9-19
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 %
%%%%%%%
Pablo Escobar
Pablo Escobar 2016-9-19
PFA Answers

类别

在 Help Center 和 File Exchange 中查找有关 General Applications 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by