proc3.m

Matlab学习课件

% Script file to run a simulation of the three-tank process system with 
%moveable filling carriage
% This version uses the single-thread, sequential scheduler
% Created June, 1995, D. M. Auslander
% Modified to match new template: 8/2/95, DMA
% File: proc3.m

procglbl	% Read in global definitions

% System setup and initialization

ncall = 0;	% Used to check how many evaluations the ODE integrator uses
v_car = 0;
Xset_car = 0;
Xnext_car = 0;
Acc_car = 1.5;	% Carriage move profile acceleration
Vc_car = 1.2;	% Carriage cruise velocity
pa1set = 0.3; pa2set = 0.6; pa3set = 0.4; 	% Desired proportion of component A
hlow = 3; hhigh = 4;	% Acceptable range for liquid heights in tanks
h1 = hhigh; h2 = hhigh; h3 = hhigh; 	% Initial height of liquid in tanks
ca = 1.1; cb = 0.9;	% Specific heat 
densa = 1.2; densb = 1.1;	% Density
A1 = 4; A2 = 3; A3 = 2.5;		% Tank cross-sectional area
xt1 = 1.2; xt2 = 2; xt3 = 3;		% Tank locations
temp1set = 20; temp2set = 25; temp3set = 23;
temp1 = temp1set; temp2 = temp2set; temp3 = temp3set;
temp_a = 15; temp_b = 15;
heat1 = 0; heat2 = 0; heat3 = 0;
flow_a = 0; flow_b = 0;
flow1out = 0.3; flow2out = 0.5; flow3out = 0.2;
flowmax = 3;

% Initialize the state vector
% State variables:
%	x(1)	carriage position, x_car
%	x(2,3,4)	mass of component A in tanks1,2,3; va1, va2 ,va3
%			(Note - these are misnamed - the v's were for volume originally!)
%	x(5,6,7)	mass of component B in tanks 1,2,3; vb1, vb2, vb3
%	x(8,9,10)	temperature in tanks 1,2,3; temp1, temp2, temp3
x = zeros(10,1);

x_car = 0;
x(1) = x_car;		% Initial value of state vector
va1 = pa1set * densa * h1 * A1; x(2) = va1;
va2 = pa2set * densa * h2 * A2; x(3) = va2;
va3 = pa3set * densa * h3 * A3; x(4) = va3;
vb1 = (1 - pa1set) * densb * h1 * A1; x(5) = vb1;
vb2 = (1 - pa2set) * densb * h2 * A2; x(6) = vb2;
vb3 = (1 - pa3set) * densb * h3 * A3; x(7) = vb3;
x(8) = temp1; x(9) = temp2; x(10) = temp3;

Kp_car = 3;
Tcar_fb = 0.2;
Tcar_fb_next = 0;	% Make sure it runs the first time
Tcar_sup = 0.2;
Tcar_sup_next = 0;
Tlevel = 0.3;
Tlevel_next = 0;

tfinal = 30;
del_t = 0.02;
del_tsim = 0.2;	% Discretization interval for simulation
ept = del_t / 100;	% Used to check for round-off errors

t_outint = 0.3;	% Time interval for outputs
		% This is particularly important for the student edition!
nouts = ceil(tfinal / t_outint);
yout = zeros(nouts + 1,4);		% Set this up fro proper output dimension
tout = zeros(nouts + 1,1);
t_outnext = 0;	% Time at which next output will be copied
iout = 1;		% Index for outputs
tstep = 0;		% Running time
ode_opts = odeset('refine',0,'reltol',1.e-3);	% Differential solver options

% Make up a task list using the string-matrix function, str2mat()
% Note that this can only take 11 strings in each call, so it must be called
% multiple times for more than 11 tasks

t1 = str2mat('car_fb','car_sup','level');
tasks = str2mat(t1);
nn = size(tasks);
ntasks = nn(1);		% Number of rows in 'tasks'

% Set initial states
Scar_fb = 0;
Scar_sup = 0;
Slevel = 0;

% Initialize all the tasks
for j = 1:ntasks
	tsk = tasks(j,:);	% name of j-th task
   if j==1 
      feval('car_fb'); 
   end;
   if j==2 
      feval('car_sup'); 
   end;
   if j==3 
      feval('level'); 
   end;
   
   %feval(tsk);	% No arguments - task takes care of state
			% and when it runs internally.
end

i_audit = 1;	% Index for transition audit trail
tbegin = get_time;		% Time at the beginning of simulation step
 % Step out the solution
while tstep <= tfinal
	% Compute the tank heights
	h1 = (va1 / densa + vb1 / densb) / A1;
	h2 = (va2 / densa + vb2 / densb) / A2;
	h3 = (va3 / densa + vb3 / densb) / A3;

	if (tstep + ept) >= t_outnext
		t_outnext = t_outnext + t_outint;	% Time for next output
		yout(iout,1) = x(1);
		yout(iout,2) = h1;
		yout(iout,3) = h2;
		yout(iout,4) = h3;
		tout(iout) = tstep;
		iout = iout + 1;
	end
	
	% Run the control tasks

	for j = 1:ntasks
		tsk = tasks(j,:);	% name of j-th task
      if j==1 
  	   	 feval('car_fb'); 
   	end;
   	if j==2 
      	feval('car_sup'); 
   	end;
   	if j==3 
      	feval('level'); 
   	end;

      %		feval(tsk);	% No arguments - task takes care of state
				% and when it runs internally.
		inc_time;	% Increment time
		del = get_time - tbegin; % Time since last simulation
		if del >= del_tsim	% If true, do a simulation step
			tbegin = get_time;
			% This shows alternate ODE solvers -- odeeul is 
			% lightning fast but has absolutely no error control!
			% ode23 uses the new Odesuite, modified to allow 
			% 'hfirst' input. Comment out the unused one.
			%ode_opts = odeset(ode_opts,'hfirst',del);
				%Specify the intial step size to avoid
				% thrashing in the ode solver
			%[tt,xx,stats] = ode23('proc_sys',del,x,ode_opts); 
			[tt,xx,stats] = odeeul('proc_sys',del,x,ode_opts); 
			% simulate control object over one discrete time step
			x = xx;	% Update state vector 	
			x_car = x(1);
		end
	tstep = get_time; % Update time
	end
end
yout(iout,1) = x(1);
yout(iout,2) = h1;
yout(iout,3) = h2;
yout(iout,4) = h3;
tout(iout) = tstep;
% Plot results
plot(tout,yout(:,1),'k',tout,yout(:,2),'k',tout,yout(:,3),'k--',tout,yout(:,4),'k-.');
xlabel('time'); title('Carriage Position and Liquid Levels');