User:Akohnert/Proportional-integral controller in a flow system

Consider a problem where a small tank is connected to a large pressure reservoir. The tank has a release valve to the atmosphere which is governed by a proportional-integral controller. The time behavior of this system is nearly impossible to solve analytically, but very simple to solve numerically. Given the gain (G) and reset (T) parameters of the controller, we can easily determine its output signal using the formula for a proportional integral system,

$$\mathrm{Output}=G \frac{p-p_{set}}{p_{set}}+\frac{G}{T}\int{\frac{p-p_{set}}{p_{set}}}$$

This will give us a target value for the valve position. We can interpret a signal of one to mean full open and a signal of zero to mean fully closed. We will assume a linear dependence between the valve position and flow coefficient. The valve coefficient ($$c_v$$) is used in the Bernoulli equation to determine the flow rate out of the tank.

$$Q=c_v\sqrt{\frac{P_1-P_2}{\gamma}}$$

We have to be careful to limit the valve coefficient to 1 (fully open) even if the output signal rises above 1. The same equation can be used with the known inlet valve coefficient to determine the inlet flow from the pressure reservoir. Using the two flow rates, we can find the change in tank pressure using a "tank capacitance factor". This is a constant that relates the difference in inlet and outlet flows to the differential pressure change in the tank,

$$\frac{\delta P}{\delta t} = C \left( Q_i-Q_o \right)$$

This is all the requisite information to begin writing the program. We need to know the pressure as time progresses, so the program will move through small increments in time. The structure will be as follows:
 * Calculate the flow rates based on existing value of tank pressure
 * Calculate the change in tank pressure due to flow rate differences
 * Using the new tank pressure, calculate the output signal
 * Adjust the valve position to match the output signal
 * Increment the time
 * Tabulate time, valve position, and tank pressure
 * Repeat

An example of this code is provided below for Mathematica. The first block sets the initial conditions, namely a reservoir pressure of 100 psi, a target pressure of 80 psi, a tank capacitance of 0.5 psi per gpm per second,and a time step of 0.05 seconds. cvo = 0; pi = 100; po = 0; psp = 80; G = 2; T = 1; $$\Delta$$ t = .05; $$\epsilon$$ = {}; Op = 0; cvi = 1.484; P = {100}; t = {0}; Qi = 0; Qo = 0; c = .5; $$\gamma$$ = 62.4; currentP= 0; cvolist = {0};

The second block represents the steps in the loop outlined above. Notice that the rate of change in the valve position has been capped to reflect real limits on the speed of opening and closing valves. For[i = 1, i <= 200$$\Delta$$t, i = i + 1, {  (*Select most recent pressure*) currentP = P[[i] ];  (*Determine the change in pressure tank capacitance factor and current flow rates*)   dp=$$\Delta$$t/c*(cvi*Sqrt[(pi-currentP)/$$\gamma$$]-cvo*Sqrt[(currentP-po)/$$\gamma$$]);   (*Tabulate values for error, new pressure, and time*)   AppendTo[$$\epsilon$$, (currentP - psp)/psp];   AppendTo[P, currentP + dp];   AppendTo[t, $$\Delta$$t*(i)];   (*Determine the controller output signal using the given formula*)      Op = G*$$\epsilon$$[[i] ] + G $$\Delta$$t/T*Total[$$\epsilon$$];   (*Adjust valve position in accordance with the output signal*)   If[Op > cvo + .01, cvo = cvo + .25 $$\Delta$$t];   If[Op < cvo - .01, cvo = cvo - .25 $$\Delta$$t];   (*Ensure real limits on valve position*)   If[cvo > 1, cvo = 1];   If[cvo < 0, cvo = 0];   (*tablulate valve position*)   AppendTo[cvolist, cvo];  };]