From simulation to computer-aided design of control systems

Control systems are among the essential engineering achievements. We consider them completely natural, but control system designers have to struggle before they succeed in their job. Yes, there are available many simulation systems, like Simulink (under Matlab). If the programmer needs to implement more than just simulation, for example, an application to automatically analyze a controlled system and design an optimal controller, using a suitable general-purpose programming language, is more desirable. Most applications run under Microsoft Windows, so let’s look how to write one in the C# language, a “native” language for such Microsoft Windows applications.

Even control system designers may find useful information here to help with simulation system programming in C#. Let’s start with the basic theory of control systems.

Control systems introduction

Figure 1 shows a basic block diagram of a generic, closed-loop control system. It consists of the controlled system or plant and the controller. The controller’s purpose is maintaining output variable at the level determined by the reference input. This is accomplished by the comparison (subtraction) of the output variable, y(t) and the reference input, r(t). The result is a regulation error, e(t), which is an input information for the controller. The controller calculates the actuating variable, u(t), which is the input to the controlled system.

While the controlled systems are predominantly analog systems (that is, they could be described by the continuous-time mathematical equations), today’s controllers are predominantly digital systems, so the block diagram of such a closed-loop control system can be represented by Figure 2. Two new components are the digital to analog (D/A) and analog to digital (A/D) converters.  Only two signals, u(t) and y(t) are of the continuous-time character. The remaining signals are in the discrete form. So, the entire signal processing starting with the comparison of r[n] and y[n] can be accomplished by a suitable digital system (a microcontroller, for example).

Control system design: System analysis

The first phase of a control system design is the controlled system (plant) analysis. You need to extract as much information about the controlled system as possible. For the control purposes, area of interest is the transfer function of the system, preferably in a well-known (s-domain) form:

If you are designing the controlled system yourself, then you might be able to model (describe) it with a set of differential equations. It sounds easy, but even for the simplest controlled system, like a servomechanism with a dc servo motor, it can take enormous amount of time to finalize its mathematical description.

If you cannot mathematically describe your controlled system, treat it as a “black box” and find its transfer function by a different way. The best way is to apply a step function at the input to the system and to capture and plot its time response.  What you will get will be very likely a “S” shaped curve as it is shown in Figure 3.

Now plot a tangent line going through the inflection point of the curve as it is shown in Figure 3. On the time (t) axis you should be able to mark two distinguished line segments, T_u and T_n. From that ratio, users will be able to determine the order of the transfer function, and, in the case of the most common, second-order transfer function, you can find out both time constants. Controlled system analysis stops here, though you might find that your controlled system is equivalent to the second-order system with the following transfer function (where the gain, K, is one for simplicity).

Regarding the gain value, we need to look at it as to a relative value. System has the gain of 1 if, for example, 50% of the input signal/variable causes the system output to achieve exactly 50%. And this is how a well-designed control system should behave. If the system output is lower than those 50%, the system is under-sized (K  < 1), it will never achieve its maximum output. The opposite case is an over-sized system (K > 1), which will reach the full output even if the input is below its maximum.

So, how can we write a computer program for the simulation of controlled systems? Well, at first we need to find the way how to describe such continuous-time systems with the discrete-time tools. Any such system can be described as a recursive filter, or Infinite Impulse Response (IIR) filter. Output of such a filter can be mathematically described as the sum of the following difference equation:

y[n] = a₀x[n] + a₁x[n-1] + a₂x[n-2] + … + b₁y[n-1] + b₂y[n-2] + …

where y[n], y[n-1],… are the system (filter) outputs at the current timestamp, n, previous timestamp, n-1 etc., and  x[n],  x[n-1], …  are the input values to the system at the current timestamp, n, previous timestamp, n-1 etc. a₀, a₁,… b₂,… are the weight coefficients forming the convolution kernel. They are usually referred to as the recursion coefficients.

Output of a single pole (that is, first-order) low-pass filter is just:

y[n] = a₀x[n] + b₁y[n-1]

For the synthesis of control systems is useful to know the transfer function (that is, y/x  ratio), which in the case of this simple, first-order low-pass filter, is:

H[z] = a₀ / (1 – b₁z^-1)                            (expressed in the Z transform/domain)

Now we need to find the relation between those two coefficients, a₀ and b1, and the time constant of a corresponding single pole filter. In an RC (analog), filter the R*C product represents the time constant of a filter, that is, the time (in seconds), which takes an RC filter output to rise to 63.2% of its input (step function) value, or decay to 36.8% of its original step function value (if such a RC filter is a high-pass filter). This rise or decay follows a simple exponential function. For a single pole low-pass filter there is a simple rule for selection of those two coefficients:

a₀ = 1 – x
b1 = x

…where x is a value between 0 and 1, and it actually represents the amount of decay (or rise) between two samples. The higher number, the slower decay/rise. In a discrete-time representation we are rather interested in the representation of this decay/rise in terms of samples. For example, we want to create a recursive filter with the “time constant” of 10 samples, which means that after exactly ten samples (steps or iterations) the filter output will decay/rise exactly by those 63.2%. Here is a simple formula, which expresses the decay coefficient, x, by a number of samples, d :

x= e^-1/d

Now, we can write a code representing such a recursive, low-pass filter. While you can use any programming language, the best programming language for such purposes is any modern, object-oriented language. But for the reason mentioned earlier let’s look at the code written in C#:

When we instantiate an object of the FirstOrderSystem class, we will prepare for it the a₀ and b1 coefficients based on the desired number of samples, d, representing the time constant of the filter. When we then call this object’s Response() method with an argument representing the step function value, it will return the filter output after one iteration.

Users can verify how the first-order low pass filter performs by adding the following code to the Program.cs source file:

Screenshot 1 shows the console [Microsoft Windows, and others] output when you run the program above. Users can see the first-order low pass filter output after 10-th iteration is exactly 63.2 (that is, 63.2% of the settled output value).

Can we simulate, for example, the above shown second-order transfer function of the controlled system? Actually, yes, we can, we just need to connect two, first-order systems (low-pass filters) into a cascade, like this:

…where we replaced the original time constants of 4 and 16ms by 4 and 16 samples/iterations, resp. This selection makes our “sampling” rate to be exactly one sample per 1ms. It will work, just the sampling would be a little bit “rough,” as to have four samples of the time constant is too few samples. Ten times higher sampling rate would be better. But then you wouldn’t want to print a thousand values in the output console and analyze them, you would prefer to see the output values in a graphical form. So, instead of a console application, users need to write the program as a Microsoft Windows Forms application and to try some graphical representation of the results.

Before going any further, understand how many samples (iteration steps) are needed to take from the analyzed systems. For such exponentially behaving systems, it takes around 3 to 5 times more time than is their (total) time constant, to exhibit their stabilized, or settled output. Our second-order system has a combined/total time constant 20ms (4 + 16), so capturing a 100 ms interval is what we need. And taking 10 samples from 1ms interval is reasonable as well. So, count on taking around 1000 samples.

While we can simulate any second-order system by a cascade of two, first-order systems, it is not a very elegant way, as we can easily write a code of a second-order system – low-pass filter, as for example the following:

This class constructor gets two parameters, d1 and d2 when the object is instantiated, and inside the a₀ and b1 coefficients of two such first-order systems are prepared and calculated, named as a and b for the first system, and c and d for the second system. Then their combined coefficients are created and calculated to implement the second-order transfer function:

H[z] = a₀ / (1 – b₁z^-1) * c₀ / (1 – d₁z^-1) = a₀ c₀ / (1 – (b₁+d₁)z^-1  + b₁d₁z^-2  )

Control system synthesis: From s-domain to z-domain

The controlled process is often a continuous process; all the variables are continuous functions of time. On the other side the digital closed-loop controllers, instead of using continuous-time signals (reference inputs, actual values, etc.) they process only digitized samples of those analog signals.  Such systems are best described in the z-domain. However, the entire control system synthesis can be done in the s-domain, and switching to the z-domain can be done in the implementation phase.

From the control theory you can remember proportional, integral, derivative (PID) compensation is one of the most common forms of the closed loop control. Why it is so popular? In most applications, the controlled process can be expressed by a first or a second-order transfer function. The PID controller can cancel, or at least significantly compensate exactly two poles of the transfer function.

Figure 4 shows such a continuous-time (analog) PID controller. The controller calculates the regulation error, e(t), as a difference between the reference value or set point, r(t), and the actual output value, y(t). Then the proportional component of the controller multiplies this error signal by a proportional constant, K_P. The integration component calculates an integral of e(t) and multiplies it by an integral constant, K_I. The derivative member calculates a derivative of e(t) and multiplies it by a derivative constant, K_D. The sum of these three components creates the actuating value, u(t), which is applied to the controlled system.

This is the differential equation of the PID control block output as a function of the regulation error:

The transfer function of the PID controller,   , in the s-domain is expressed as

Now, let’s have a controlled system, which can be approximated by a second-order transfer function

Then the open loop transfer function of the entire control system,  , equals to

And here you can see, why such a PID controller is so important. You can optimally tune its parameters, like the  ratio can match the value of τ₁ + τ₂ sum and the  ratio can match the value of τ₁ * τ₂ product, so the open loop transfer function will be reduced just to:as both poles of the transfer function have been “cancelled.” In such an optimally tuned control system the closed loop transfer function becomes just a first-order, low-pass filter system with the time constant, τ, equivalent to (if we assume the controlled system gain, K, to be 1). This is very good, because we can achieve much faster time response from the closed loop control system (if you make Κ₁ adequately large) than what is the response of the controlled system itself.

Control system implementation, simulation

The analog PID controller can be modified to be implemented as a discrete-time control system, as it is not difficult to rewrite the differential equation of the PID controller into its difference form:

where u[n] is the actuating value at the present time n, e[n] is the regulation error at the time n and

e[n-1] is the regulation error at the previous sample time, n-1. T is the time period of the sampling. The same time period, T, is used for the processing, such as for the u[n] calculation.

For the practical applications certain modifications of this equation are needed. At first, the integral member will better use the trapezoidal approximation of integration instead of the rectangular one. Then the derivative member tends to be noisy. To reduce noise, we can use more than two (for example four) consecutive samples of the regulation error and put them through a tiny (4-tap) FIR filter. Here is the modified PID formula:

The difference equation above can be easily implemented in C#, and the PIDControlBlock class can look like this, for example:

As we already have the controlled system (plant) class designed, we can put together those two objects and create the final, ClosedLoopSystem class:

Please, notice that the above shown ClosedLoopSystem class is ready to be used not only with the SecondOrderSystem objects, but with another, commonly used, FirstOrderAstatic systems/objects as well. They represent entire class of the servomechanisms. You can create the FirstOrderAstatic class combining a FirstOrderSystem class code with an integrator. Here is its transfer function:

And here is the code in the Main() function, with which users can simulate the entire closed-loop control system. The plant object is represented as a second-order system with the time constants 250ms and 100ms (or iteration/samples). The PID object is instantiated with the optimal, Kp, Ki, Kd, coefficients, and with the Ki value selected such, that it will make the closed-loop response 5-times faster than the combined time constant of the plant.

Now run the program, and this is what you will see. Actually, on the console you will not see all 1000 values, only a few hundreds of the latest values. So, to see them all, you might write them to a text file. Then you could see them all. Then try to find a value close to 63.2. You will find it at a line 70 or so. What does it mean? That our second-order controlled system with the combined time constant 350 (250 + 100), now behaves as a system with the time constant around 70. Indeed, our closed-loop control system has the (total) time constant five times smaller.

Computer-aided design application

Once you understand the nature of control systems analysis and synthesis, and you know how to represent/simulate it on a computer, it is not so difficult to write a more complex computer program, one which would serve as a computer-aided, control system design application.

Automating the control system analysis is not so difficult. At first, you need to get into your program the data representing the open loop transfer function, such as the “S” curve values. They can be saved in a text file. Then calculating the derivative of the S-curve data will help you to calculate the T_u and T_n segments.

Then by following a suitable look-up table and interpolation, your program can calculate the time constants of the analyzed controlled system. Author wrote such a computer-aided control systems design application, and here you can see how precisely those time constants could be calculated. Screen shot 2 shows the analysis result of a second-order system having two 250 ms time constants:

Screenshot 3: First-order astatic system analysis result is shown in a perfect matching. Courtesy: Peter Galan, a retired control software engineer[/caption]

Once the application can load and analyze data representing the transfer function of a controlled system, users can start designing the “control system synthesis” part of the program.

In that part, users will leverage their knowledge of tuning the PID controller, which will be the base of the entire control system synthesis. But that’s not enough. While you can perfectly tune any PID controller, such a control system has a limited practical usage. Why? Because in real systems your real PID controller won’t be able to provide (almost) unlimited actuating variable, u(t), which the optimal control requires. Yes, an optimal controller will “generate” even 100 times higher values, than what the real controlled system can accept. Imagine a temperature control system using a 10kW heater, which would have applied, though for a short time, a power corresponding to 1MW! So, a “realistic” PID controller has to take into account (and limit!) the actuating variable, u(t) to some maximum value. Author’s application assumes the provided transfer function (“S” curve) was taken when the input variable (step function) was at 50% of its maximum value. At the gain of 1 such system output will reach 50% of its maximum. So, what is a reasonable limit of the actuating variable – output of the PID controller? A very reasonable selection is to use 2-times of the last (stabilized) value found in the captured “S” curve data.

In Screenshot 4, users can see the transfer functions of an ideal (green) and realistic (red) closed-loop control system. Limited actuating variable causes a huge overshot, which is usually unacceptable, so some other improvements are needed. One solution is to add a Feed-forward member to the control system as it is shown in the following screen shot as well. And the result is a transfer function shown in blue. Not ideal, but still much better than the open-loop response of the system (black).

Screenshot 5: First-order astatic control system synthesis is show from a time-optimal (bang-bang) controller. Courtesy: Peter Galan, a retired control software engineer[/caption]

What else can be added to such a computer-aided design application? Look at the following Screenshot 6. It shows the “Manual Tuning” form, in which the user can change the values of the control parameters, and concurrently their immediate impact is seen. This can be very useful, as the user can realize much better how the individual control parameters affect the entire control process. And with a fine tuning the user can obtain even better more “likeable” results.

Control systems can be simulated in C# or Python

Complex systems like the control systems can be relatively easily simulated by a suitable programming language like C#. If you are not a professional programmer, maybe Python (which can be easily learned) may be a better choice. Simulation is the simpler part of such a project. System analysis is a little more demanding. You would need to analyze the “S” curve data, and with the help the look-up-table and interpolation, the time constants can be found. More demanding still might be the results representation. Fortunately, Microsoft Windows (with its .NET platform) offers a rich library of the supporting classes and other programming items for such data visualization. Here is what you would need to add to your program:

using System.Windows.Forms.DataVisualization.Charting;

Then the task becomes easy. On the internet, people can find dozens of examples how to use such data visualization. Good luck in future endeavors in creating computer-aided, control systems design applications.

Peter Galan is a retired control software engineer. Edited by Mark T. Hoske, content manager, Control Engineering, CFE Media, mhoske@cfemedia.com.

KEYWORDS: C# control system programming, programming simulation, PID

CONSIDER THIS

Are you using the right tools for programming a control system design?

Do you have experience and expertise with the topics mentioned in this content? You should consider contributing to our CFE Media editorial team and getting the recognition you and your company deserve. Click here to start this process.