I am simulating the height of the reservoir in a hydro power plant using Matlab, without Simulink. The radial gates of the power plant are used to control the height of the water in the reservoir. This is done by implementing a PI controller on the radial gates. First, I am trying to let things work with only one radial gate.
The flow into the reservoir is chosen to be
$$Q_{in} = 400 \,\,m^3/s. \qquad \qquad (1)$$
The output flow through a radial gate is given by
$$Q_{out} = \frac{dV}{dt} = C_g(D) \cdot D \cdot L \cdot \sqrt{2\cdot g \cdot H_g} \qquad \qquad (2)$$
where
$$H_g = H - \frac{D}{2} \qquad \qquad (3)$$.
Rewriting eq. (2) with respect to (3) gives $$Q_{out} = C_g(D) \cdot D \cdot L \cdot \sqrt{2g (H - \frac{D}{2})} \\ = C_g(D) \cdot D \cdot L \cdot \sqrt{g} \cdot \sqrt{2H-D}. \qquad \qquad (4)$$ Grouping all the constants as $C_1 = L \cdot \sqrt{g}$ gives
$$Q_{out} = C_g(D) \cdot D \cdot C_1 \cdot \sqrt{2H-D} \qquad \qquad (5)$$
Eq. (5) includes following variables and constants:
- $C_g(D)$ is a function dependent on the gate opening $D$.
- $D$ is the gate opening. It's interval is 0 - 8.7 meters.
- $C_1$ consists of $L$ which is the width of the gate and gravity $g$.
- $H$ is the height of the reservoir, calculated from the edge of the lower part of the gate where $D = 0$.
Furthermore, I am using a PI controller to control the height $H$. I am using a for-loop in Matlab to simulate this problem. The discrete PI controller I use is calculated as follows.
The error in Matlab is calculated by $$error(1,i) = H(1,i) - H_d(1,i)$$ where $H$ is the height and $H_d$ is the reference value or the set point. The height of the reservoir should always be at $116 m$ above sea level. In the program, I subtract the lower edge of the radial gate, where $D = 0m$, which is $H_{edge} = 110m$, so the radial gate is placed very high above sea level. I also choose the initial height in the simulation to be $H_{initial} = 115.8 m$, therefore, $$H(1,1) = H_{initial} - H_{edge} = 115.8 - 110 = 5.8m \\ H_d(1,1) = H_{d} - H_{edge} = 116 - 110 = 6m.$$ This gives an initial error of $error = -0.2m$.
Next thing I do is to calculate the $P$ and the $I$ terms: $$PI\_P(1,i) = k_p \cdot error(1,i) \\ PI\_I(1,i) = (\frac{k_i}{T_i}) \cdot (PI\_I(1,i-1) + 0.5 \cdot (error(1,i) + error(1,i-1))\cdot dt) - intWindup(1,i)$$
The integral term is actually calculated in the program as
if i == 1
PI_I(1,i) = ((ki/Ti)*error(1,i)*dt) - intWindup(1,i);
else
PI_I(1,i) = ((ki/Ti)*(PI_I(1,i-1) + 0.5*(error(1,i) + error(1,i-1))*dt)) - intWindup(1,i);
end
The controller output is then calculated as the sum of the two terms: $$ PI\_u(1,i) = PI\_P(1,i) + PI\_I(1,i) $$
Since the radial gate can only vary between 0m up to 8.7m, I used a integral windup method called "backwards calculation" so that the integrator will not "windup". This is shown in the image below:
Next thing I do is calculate the saturation $$u\_sat(1,i) = CheckSaturation(PI\_u(1,i),Dmax,Dmin)$$ The function I wrote is given below
function u_sat = CheckSaturation(u,gateMax,gateMin)
% Controller saturation check
if u >= gateMax
u_sat = gateMax;
elseif u <= gateMin
u_sat = gateMin;
else
u_sat = u;
end
end
After this, I calculate the integral windup as $$intWindup(1,i+1) = intWindupBackCalculation(PI\_u(1,i),u\_sat(1,i),k\_lim)$$ where $k\_lim = 1$. The function I wrote is given by
function val = intWindupBackCalculation(u,u_sat,k_lim)
diff = u - u_sat;
val = diff * k_lim;
end
Now finally, I update the model of the reservoir with the following code
Cg1(1,i) = Cg_D(D1(1,i),Dmin,Dmax);
% The physical system
Qout_gate1(1,i) = Cg1(1,i) * C1 * D1(1,i) * sqrt(2*H(1,i) - D1(1,i));
Qout(1,i) = Qout_gate1(1,i) + Qout_m;
H_dot(1,i) = (1/A(1,i)) * (Qin(1,i) - Qout(1,i));
H(1,i+1) = H(1,i) + 0.5 * (H_dot(1,i) + H_dot(1,i+1)) * dt;
This is simply calculation of eq. (2) from above, $Qout\_m = 0$ is the flow through the machines which is excluded in the calculations. The differential equation used to calculate the change of height is $$ \dot H = \frac{1}{A} (Q_{in} - Q_{out}) $$ then in the last line of the code, I use discrete integral to calculate the height.
My question starts now.
Below is a plot where I choose $k_p = 40$ and $k_i = 0.1$
- $H$ is the change of height per second
- $H_d$ is the desired height, or the reference level, 116m
- $D1$ is the radial gate opening in meters
- error is simply the error, the blue line, while the integral windup is the red line. it is also possible to see the proportional and integral gains in the legend.
- Qin is the red line and Qout is the blue line
- $PI\_P$ is the proportional term of the controller
- $PI\_I$ is the integral term of the controller
In these figures above, it can be seen that the controller tries to keep the water height at 116m but there is this settling error, I think it's related to slowness of the integral term. It can also be seen that the radial gates opens up and it tries to regulate the water heigt. The error can be seen as very small related to the integral windup, but the error is there and it's not always zero. The flow in vs. flow out starts to be equal slowly. On the last figure, it can be clearly seen that the proportional term and the integral term is working.
Now I choose $k_p = 40$ again, but the integral gain $k_i = 1$
In the figures above, the settling error is much smaller now and the response is much faster. The integrator is doing much more now, but not enough. I want the water level to be at 116m above sea level with only +/- 5cm. In real life, this system is very slow. The reservoir is huge, its area is $A = 3.8 \cdot 1000^2 \, \, m^2$ or $3.8 \,\, km^2$.
Now I want to increase the integral gain, and go above 1. So I choose $k_p = 40$ and $k_i = 2$.
In the figures above, it can clearly be seen that everything is very slow now which is not acceptable.
Now I increase the integral term to $k_i = 5$ and everything goes unstable and crashes as you can see in the figures below.
My question is: What can I do about this problem? Do you have any idea of what I'm doing wrong or what I could do better? The integrator is being a little bit of problem for me and I thing it's a bit strange that the system crashes.
The whole code is below:
t0 = 0.0; % [s] Initial time
days = 0;
hours = 8;
minutes = 60;
seconds = 60;
F_freq = 50; % Frequency (increase -> smaller dt, decrease greater dt)
Fs = 2*F_freq; % Samples per second (Sample frequency Fs must be at least two times the frequency F_freq to avoid aliasing)
dt = 1/Fs; % Seconds per sample
if days == 0 && hours ~= 0 && minutes ~= 0 && seconds ~= 0
t_f = seconds*minutes*hours; % [s] Final time
elseif days == 0 && hours == 0 && minutes ~= 0 && seconds ~= 0
t_f = seconds*minutes; % [s] Final time
elseif days == 0 && hours == 0 && minutes == 0 && seconds ~= 0
t_f = seconds; % [s] Final time
else
t_f = seconds*minutes*hours*days; % [s] Final time
end
t = (t0:dt:t_f-dt)'; % [s] Sample instants
T_plot = t';
Nsim = length(t);
% A = 3.8 * 1000^2; % [m^2] Area of Hagalón, "equivalent to 3.8 km^2"
% Qout_machine = 352; % [m^3/s] Max flow through machines.
% Qin(1,1:Nsim) = [200*ones(1,round(Nsim/8)) 300*ones(1,round(Nsim/8)) 400*ones(1,round(Nsim/8)) 500*ones(1,round(Nsim/8))...
% 600*ones(1,round(Nsim/8)) 700*ones(1,round(Nsim/8)) 600*ones(1,round(Nsim/8)) 500*ones(1,round(Nsim/8))]; % [m^3/s] Flow into Hagalón
Qin(1,1:Nsim) = 400;
H_edge = 110;
H_initial = 115.8 - H_edge; % [m] Initial height of the lagoon
% Hd = [116*ones(1,round(Nsim/2)) 116.4*ones(1,round(Nsim/2))] - H_edge;
% [m] H desired (reference value)
Hd(1,1:Nsim) = 116 - H_edge; % [m] H desired (reference value)
% Parameters for radial gate
L = 12; % [m]
g_sqrt = sqrt(9.82); % [m/s^2]
C1 = L*g_sqrt; % Group C1 as a constant [m^2/s^2]
Dmax = 8.7; % [m] Maximum opening D of a radial gate
% Dmax = 15; % [m] Maximum opening D of a radial gate
Dmin = 0; % [m] Minimum opening D of a radial gate
a_slope = 8.7/522;
% Model
H_dot(1,1:Nsim) = 0;
H(1,1:Nsim) = H_initial;
D1(1,1:Nsim) = 0;
Cg1(1,1:Nsim) = Cg_D(D1(1,1),Dmin,Dmax);
A(1,1:Nsim) = 3.8 * 1000^2;
Qout_m = 0; % Flow out of two machines. Max flow is Qout_m = 352 m^3/s
Qout_gate1(1,1:Nsim)= 0;
Qout_gate1(1,1)= Cg1(1,1) * C1 * D1(1,1) * sqrt(2*H(1,1) - D1(1,1));
Qout(1,1:Nsim) = 0;
% PI
error(1,1:Nsim) = H(1,1) - Hd(1,1);
PI_P(1,1:Nsim) = 0;
% PI_I(1,1:Nsim) = 0.000518663;
PI_I(1,1:Nsim) = 0;
PI_u(1,1:Nsim) = 0;
u_sat(1,1:Nsim) = 0;
% P controller for the radial gate opening
error_gate(1,1:Nsim) = 0;
PI_P_gate(1,1:Nsim) = 0;
PI_I_gate(1,1:Nsim) = 0;
PI_u_gate(1,1:Nsim) = 0;
u_sat_gate(1,1:Nsim) = 0;
% P and I gain coefficentes for the height of the reservoir
kp = 40; % Proportional Gain
% tau_i = seconds*minutes; % Time constant for I term
% ki = 1/46104;
ki = 5;
k_lim = 1;
% ki = 1/tau_i; % Integral Gain
% ki=0;
% ki=kp;
% Ti = 1/dt; % Time constant for integral term
Ti = 1;
% P gain coefficient for the gate opening
kp_gate = 1;
% Anti-Reset Windup switch coefficient
kc(1,1:Nsim) = 1;
intWindup(1,1:Nsim) = 0;
kc_gate(1,1:Nsim) = 1;
for i = 1 : Nsim-1
% PI controller for the height in the reservoir
error(1,i) = H(1,i) - Hd(1,i); % Correct one
PI_P(1,i) = kp*error(1,i);
if i == 1
PI_I(1,i) = ((ki/Ti)*error(1,i)*dt) - intWindup(1,i);
else
PI_I(1,i) = ((ki/Ti)*(PI_I(1,i-1) + 0.5*(error(1,i) + error(1,i-1))*dt)) - intWindup(1,i);
end
PI_u(1,i) = PI_P(1,i) + PI_I(1,i);
% Check saturation of the radial gate
u_sat(1,i) = CheckSaturation(PI_u(1,i),Dmax,Dmin);
% Calculage anti-reset windup calculation using clamping
% kc(1,i+1) = intWindupClamping(error(1,i+1),PI_u(1,i),u_sat(1,i));
intWindup(1,i+1) = intWindupBackCalculation(PI_u(1,i),u_sat(1,i),k_lim);
% Feedback for the opening of the radial gate (slowness)
if i == 1
error_gate(1,i) = 0;
D1(1,i) = kp_gate*a_slope*error_gate(1,i);
else
error_gate(1,i) = u_sat(1,i) - D1(1,i-1);
D1(1,i) = D1(1,i-1) + kp_gate*a_slope*error_gate(1,i);
end
% D1(1,i) = u_sat(1,i);
Cg1(1,i) = Cg_D(D1(1,i),Dmin,Dmax);
% Cg1(1,i) = 0.7;
% The physical system
Qout_gate1(1,i) = Cg1(1,i) * C1 * D1(1,i) * sqrt(2*H(1,i) - D1(1,i));
Qout(1,i) = Qout_gate1(1,i) + Qout_m;
H_dot(1,i) = (1/A(1,i)) * (Qin(1,i) - Qout(1,i));
H(1,i+1) = H(1,i) + 0.5 * (H_dot(1,i) + H_dot(1,i+1)) * dt;
% A(1,i) = AreaReservoir(H(1,i+1));
end
PI_P(1,i+1) = kp*error(1,i+1);
Qout_gate1(1,i+1) = Cg1(1,i) * C1 * D1(1,i+1) * sqrt(2*H(1,i) - D1(1,i+1));
figure; % Change of height
%sgtitle('Change of height per sec - Tank')
sgtitle('Change of height per sec - Tank')
subplot(4,1,1)
plot(T_plot,H + H_edge,"LineWidth",2); hold on;
plot(T_plot,Hd + H_edge,"--","LineWidth",2);
legend("H","Hd - reference");
ylabel('H [m]'); xlabel('Time [s]')
xlim([t0 t_f]); grid;
% ylim([0 120]); grid;
subplot(4,1,2)
plot(T_plot,D1,"LineWidth",2); hold on;
legend('D1')
ylabel('D1 [m]'); xlabel('Time [s]')
xlim([t0 t_f]); grid;
subplot(4,1,3)
plot(T_plot,error,"LineWidth",2); hold on;
% plot(T_plot,kc,"LineWidth",2);
plot(T_plot,intWindup,"LineWidth",2);
ylabel('error [m]'); xlabel('Time [s]')
xlim([t0 t_f]); grid;
comment_1 = sprintf("error - kp = %.6f",kp);
comment_2 = sprintf("Windup - ki = %.6f",ki);
legend(comment_1,comment_2)
subplot(4,1,4)
plot(T_plot,Qout,"LineWidth",2); hold on;
plot(T_plot,Qin,"LineWidth",2);
ylabel('Q [m]'); xlabel('Time [s] (2 dagar)')
xlim([t0 t_f]); ylim([0 (max(Qin)+100)]); grid;
legend("Qout [m^3/s] & D = 5m", "Qin [m^3/s] & D = 5m")
% title("Qin & Qout vs. time")
figure; % Change of height
%sgtitle('Change of height per sec - Tank')
sgtitle('PI controller - Tank')
subplot(2,1,1)
plot(T_plot,PI_P,"LineWidth",2);
ylabel('PI_P [m]'); xlabel('Time [s]')
xlim([t0 t_f]); grid;
subplot(2,1,2)
plot(T_plot,PI_I,"LineWidth",2);
ylabel('PI_I [m]'); xlabel('Time [s]')
xlim([t0 t_f]); grid;
function Cg_vs_D = Cg_D(D,Dmin,Dmax)
% Calculation of Cg variable of a ratial gate
n = 7; % Degree of the polynomial, got it from the interpolation process for Cg vs. D
if D < Dmin
Cg_vs_D = 0.7743;
elseif D > Dmax
Cg_vs_D = 0.6513;
else
Cg_vs_D = -5.1854385647666910937633559519621684330559219233692e-07*D^n ...
+ 2.3009388243133869431924706794312385227385675534606e-05*D^(n-1) ...
- 0.00039326051774348001841691280233703764679376035928726*D^(n-2) ...
+ 0.0032245124407480371502010552120509601081721484661102*D^(n-3) ...
- 0.012403237757656533982175695030036877142265439033508*D^(n-4) ...
+ 0.01661740981151977464280733443047211039811372756958*D^(n-5) ...
- 0.012264797555097968484449921788836945779621601104736*D^(n-6) ...
+ 0.77434799435072565465532079542754217982292175292969;
end
end
function u_sat = CheckSaturation(u,gateMax,gateMin)
% Controller saturation check
if u >= gateMax
u_sat = gateMax;
elseif u <= gateMin
u_sat = gateMin;
else
u_sat = u;
end
end
function val = intWindupBackCalculation(u,u_sat,k_lim)
diff = u - u_sat;
val = diff * k_lim;
end