I have been playing around with level control tuning optimisation using Scilab. What I have come up with is a draft Scilab program that will find the optimum PID settings to minimise a cost function with configurable weighting coefficients so that you can specify the relative importance of: a) Max level error on inflow disturbance b) Sum of error c) Small valve movements all on an inflow disturbance. (I should really add setpoint responses too)
Please note that I haven't properly tested it, there are most probably bugs in.
The using of integral gain instead if time constant is because you get better optimiser performance like that.
To use it, copy from below the "level Optimse.sce" and "level.sci" into two separate files and change the Scilab directory to there and then "exec" the "level Optimse.sce" from the menu in Scilab. It will then prompt you for values and return the optimised tuning coefficients.
Regards Pieter Steenekamp
*************************************************
- Fiirst the "level Optimise.sce" program *
************************************************* // This program finds the PID settings to optimise the response // on a level controller with significant dead time in the actuator
// it uses the SciLab functions "NDcost" and "optim" // to do the optimisation
clear; plotON=0; exec('level.sci') Thold=evstr(x_dialog('Tank holdup in seconds ?','50')); Tpdt =evstr(x_dialog('actuator dead time ?','10')); maxErrorWeight=evstr(x_dialog('maxErrorWeight ?','10')); sumAbsErrorWeight=evstr(x_dialog('sumAbsErrorWeight ?','1')); sumDeltaCoSquaredWeight=evstr(x_dialog('sumDeltaCoSquaredWeight ?','1'));
// initial tuning values: x(1)= -1; // for proportional gain x(2)= 0; // for integral gain x(3)= 0; // for derivative gain
printf(' Maybe now is a good time to have coffee, because this optim- function is sometimes slow \n'); [f,xopt,gopt]=optim(list(NDcost,level),x)
Kc=xopt(1) Ki=xopt(2) Kd=xopt(3)
plotON=1; g=level(xopt);
*************************************************
- Now for the "level.sci" function *
*************************************************
function f = level(x)
// The basic intent of this function is to be called by "optim" from outside of this to find the optimum x
// This fuction calculates the optimisation cost function as a function of Kc, Ki & Td, supplied as vector x // Thold, Tpdt, plotON,maxErrorWeight,sumAbsErrorWeight,sumDeltaCoSquaredWeight must be global variables, declared outside
// The approach on units taken is to work in percent of range and not in engineering units
// Optimisation variables: max_error = 0; sumAbsError=0; sumDeltaCoSquared=0; // To penalise fast cahnges in controller output
// Read tuning variables x Kc=x(1); Ki=x(2); Kd=x(3);
T=1; // simulation update period in seconds Kp = -T/Thold; // Tank gain - take it as negative because control on outflow PV(1)=0; // fluid level in meters CO(1)=0; // Control output in percent of full output ui(1) = 0; // integrator contribution to CO up(1) = 0; // proportional contribution to CO ud(1) = 0; // derivative contribution to CO flow_out(1)=0; // flow pumped out in cubic meters per minute flow_in = 10; new_error=0; // error in level meter old_error=0; // error in level meter
max_time = 200; // length of simulation in seconds N=round(max_time/T); // convert to sample periods
// intialise dead time buffer if Tpdt>0 then for j = 1:Tpdt, opBuffer(j) = 0; end end
for n=1:N;
t(n) = n*T; // Compute time indexes. This is needed for time index
SP(n)=0;
// calculate the new flow, as a function of controller output and dead time: if Tpdt > 0 then flow_out(n+1) = opBuffer(Tpdt) ; //pure dead time is applied to out floe else flow_out(n+1) = CO(n); end
// update dead time buffer if Tpdt>1 then for j = 0:Tpdt-2, j1 = Tpdt - j; opBuffer(j1) = opBuffer(j1-1); end end opBuffer(1) = CO(n);
// Integrate the difference between the flow in and flow out. PV(n+1) = PV(n) + Kp*(flow_in - flow_out(n))*T; if PV(n+1) > 100 then PV(n+1)=100 elseif PV(n+1) max_error then max_error = abs(new_error(n+1)); end sumAbsError = sumAbsError + abs(new_error(n+1)); up(n+1) = Kc*new_error(n+1); // calculate the proportional control output percent ui(n+1) = ui(n) + Ki*new_error(n+1); // calculate the integrator control output percent
if ui(n+1) > 100 then ui(n+1)=100 elseif ui(n+1) 100 then CO(n+1)=100 elseif CO(n+1) 0.5 then subplot(2,1,1); // Level on top plot(PV); xtitle('Level'); legend("level"); subplot(2,1,2); // Control output at bottom
plot(CO); xtitle('Control Output'); legend("CO"); end
f = maxErrorWeight*max_error + sumAbsErrorWeight*sumAbsError + sumDeltaCoSquaredWeight*sumDeltaCoSquared;
endfunction