Tuesday, December 30, 2014

Android Air Properties Calculator released

I've just uploaded the source code for the BasicAirData Air Properties Calculator. At a glance; Air Density is calculated with the CIPM 91 formula, viscosity with the Sutherland's formula. The application is Free and Open.

Figure 1. Application ScreenShoot 

Read more on the application page.

Figure 2. QR-Code and link to the Google Play Store

Saturday, December 20, 2014

High Angle of Attack Total Pressure Measurement

Measurement of airspeed at high angles of attack requires specially designed devices. In this article we introduce the Kiel total pressure probe. I'll introduce a BasicAirData design that will be tested during the 2015 year.

In a conventional Pitot tube the total pressure measurement is sensible to the angle of attack. For conventional accuracy and operating ranges, let's say for angle of attack under 10°, the impact of such a non ideal behavior is commonly acceptable. However as the angle of attack is increased beyond 20 degrees we register a remarkable loss of accuracy.

Figure 1. Kiel probe mounted on a NASA X31 Experimental plane. Source NASA Book page 7. That Kiel probe was the cause of X31 fly accident. 

To reduce probe sensitivity to the operating angle of attack and side slip a sleeve is put over a conventional total pressure tube. Have a look to the following figure
Figure 2. Kiel probe design extracted from NACA-tn-2530

The figure is a Kiel probe with a Venturi shaped shield. Nacatn-2530 reports the performances of such shield on table 1, first row. That probe can handle about 40°, appreciate the measure error vs angle of attack on page 15 Figure 5. In this Nacareport the probe is considered to operate correctly if there is a total pressure error under 1%.

Alternative designs are available, there are many industrial fields of application. For a commercial probe example have a look to figure 3.

The Flowkinetics probe is rated up to 49°, of course compared to a standard Pitot solution is remarkable; probes with similar design are used in HVAC and turbo machinery applications.

In the next period of time at BasicAirData we will test a Kiel like probe. The probe head will be as per the next figure. The inlet of the shield will have an elliptical shape.

Figure 4. BasicAirData Pitot head shield draft model

You find here below a video that includes the total pressure pipe.

Video 1. Preliminary High Angle of attack Pitot Head

Target performances requirements are, operating airspeed from 3m/s to 55m/s. Measurement error should be under 1% for angle of attack under 50°.

Follow up in this article.

Sunday, December 14, 2014

Automatic Pitot-Static Calibrator. Requirements Part 1.

If we need to assure our instrument performances it is necessary to test it. There are no shortcuts; our measurement should be compared against a reference instrument. This article is the first of a mini-series that will expose the requirements for the BasicAirData Open Pitot-Static calibrator; it is an Open project so interact freely.

A Pitot-Static probe is typically used to measure the airspeed, altitude and the climb ratio. The measurement relies on the knowledge of the total pressure and the static pressure. Total pressure is the pressure exerted by the air at the probe tip. Static pressure is the free stream pressure, referred also as barometric pressure. For further details about the Pitot-Static probe refer to the BasicAirData site.

Static pressure measurement is needed for altitude and vertical speed determination.

Total pressure measurement is necessary for airspeed calculation; often airborne instrumentation measures the differential pressure between the total pressure port and the static pressure port. In line of principle is the same of having two absolute pressure measurements, one per port. There are some practical issues with the later configuration, that aspect will not be investigated in this article.

So the first requirement for the calibrator is that it should have two pneumatic connections that can be fast fixed to the instrument under test. Both the probe and the calibrator should be accurately designed for this purpose. Commercial producers call this hardware Pitot Test Adapters and Static Test Adapters. Here below a picture from a Barfield catalog page 5

Figure 1. Barfield. Pitot Test Adapters and Static Test Adapters

This kind of adapters is not what we exactly need for the BasicAirData 8 mm probe. The 8mm probe have both the static and the total pressure ports on the tube body hence we need an adapter that can handle both the pressures at the same time. The pressure ports of the Pitot-Static Calibrator should be fitted with a removable connection hose. Same ports will be used to periodic verification of calibrator internal instrumentation.

Let's focus on the altimeter. To ensure that every physical device and software used in the altitude measurement is working properly it is necessary to check if for a known given pressure input the output is those predicted by ISA model. ISA model is a widespread used atmosphere model. The test should be carried out for a certain number of different pressures in a desired operating pressure set.

Refer to the following figure as the typical static pressure measurement equipment layout.
Figure 2. Typical static pressure measurement equipment layout

The aerodynamic behavior of the probe is not tested by a Pitot-Static calibrator. Is up to the probe designer to warrant a correct behavior, the user have the duty of carefully install the probe. By instance if the probe is not installed correctly the altitude indication during the fly can be wrong, even if the measurement equipment has been carefully tested. 
Figure 3. Typical Pitot-Static calibrator test layout

As per figure 2 the Pitot-Static calibrator should operate in closed loop. For every altitude within a test set the output of the probe is compared with the output of a reference instrument. The difference between the reference instrument output and the instrument under test output is indicated on the figure as “Error”. For each tested point the calibrator evaluates Error and if this value lies outside the required interval mark then marks as failed the test. A leading requirement is that all the test procedure should be completely automated.

Focus on the box “A”; inside this box the pressure at the static port is converted into a number. Any problem within the box “A” will be individuated by a Pitot-Static calibrator. Prior to the test all the equipment in box “A” should be checked and calibrated, otherwise the whole test will be a waste of time. In the most basic setups box “A” contains a pressure sensor connected with a pneumatic line to the static pressure port. The sensor should be calibrated to remove any treatable source of uncertainty, as can be a simple offset. I stress on the fact that if your airspeed indicator at rest on the bench top measures -3 m/s then you should take care of it prior to go further with the tests.

Given that another requirement is that Pitot-Static measurements should be available to the Pitot-Static Calibrator. With standard BAD 8mm Pitot a serial connection is available, others units may need a dedicated interface.

Now let's individuate the main requirement parameters for the Pitot-Static calibrator.
RC planes operate in the very low part of the troposphere. If you live at about 0 MSL you should be happy, check your local regulations, with an instrument calibrated within the (-50,100)m MSL range. However if you are operating from a higher altitude, for example on a hill top at 800 m MSL, then you should be able to calibrate your instrument with the same span but within the (750, 900)m MSL range. A huge part of the Earth's population lives under 4000 ASL [1 Fig, 3] 
hence it's seems adequate to set a calibration range of (-200, 4000)m MSL, using ISA model that lead to a calibration interval of (103751, 61640)Pa.
That is the pressure swing that should be managed at the static pressure port.
What are the pressures to test and the tolerances? 
FAA regulate that aspects under “FAR 43-Appendix E”

Figure 4. Tolerances table. Table I from FAR 43-Appendix E

Stated tolerance for 1000 feet (304.8m) is 30 feet (9.2m). 
We can use that tolerance values straight away, however to support GPS enhancement health and recovery systems or other high level tasks it seems more appropriate to require a tolerance of 1m (0.3048 feet). The tolerance of 1 meter is referred to the unit under test, the Pitot-Static calibrator should have and uncertainty that can handle such a requirement; an uncertainty of 0.2 m on the average of multiple measurements seems adequate. That value is the performance required to the Pitot-Static Calibrator internal static pressure reference sensor.

Using the ISA model at 0 MSL and considering an altitude variation of 1 meter we get 12 Pa of pressure variation. In the same conditions an uncertainty of 0.2m is equivalent to 2 Pa of uncertainty.

Calibration procedure should be repeated for a certain number of pressures or equivalent pressures. Best way to choose that test pressure is to fit the flight envelope of the aircraft that bear the equipment. If the aircraft is designed to operate, in respect of local regulations, between (-100,100)m MSL then is necessary to test at least the altimeter at -100m(0%), 0m(50%) and 100m(100%).If it is know that the airborne unit will operate most of the time around an altitude h then is wise to include also that altitude into the test routine.

Until now we've introduce some aspect of Pitot-Static Calibrator, with an attention to altitude calibration. Next article will cover airspeed calibration requirements.


[1] JOEL E. COHEN AND CHRISTOPHER SMALL (1998), Hypsographic demography: The distribution of human population by altitude, Proc. Natl. Acad. Sci. USA Vol. 95, pp. 14009–14014, November 1998. Retrieved online 9/12/2014.

Friday, November 21, 2014

Web Maintenance

Next month we will undertake ordinary and extraordinary Web maintenance tasks. Web disruption is not on schedule. Kindly report any anomaly. 

Take care,
BasicAirData Team

Sunday, November 2, 2014

Front removable flange test

After some "design-fix-test" iterations the new removable Pitot flange passed the pressure test; test pressure was 8 barg.

You find the free downloadable 3D models on GitHub, keep you up to date following the forum.

Have a look to the following general arrangement drawing.

Figure 1 FF-8-500 - Front Flange for 8 mm Pitot- Assembly – 20141025

Major concern was the possible leaking of the pressure lines; sealing is attained with o-rings.

Test consisted in pressurize the total pressure and the static pressure line. Pressure was applied at one line per time or on both lines at the same time. You can see the unit under test in figure 2.

Figure 2 BasicAirData Flange under test

The new flange is suitable to be used with the new air data computer, assembly view in figure 3
Flange connected by means of two pressure ports or lines; one port is for the total pressure that is routed from the probe tip and the other for the static pressure.

Figure 3 Removable Flange, Pitot and Air Data Computer, all from BasicAirData

As the flange is equipped with push-in pneumatic fittings the pressure vessel is also connected with push-in fittings; you can see the coupling detail on figure 4.

Figure 4 Push-in pneumatic fittings and coupling

Refer to figure 5.The main part of the test equipment consist in a pressure vessel. Within this vessel the pressure is adjustable by mean of a valve connected to an air compressor. On the same pressure vessel is fixed a manometer. To test a pressure line we reach the desired test pressure and then we wait 15 minutes. If after that lapse of time the manometer indication is still the same then the test for this port and pressure is considered passed.

Figure 5 Test equipment
The  maximum expected operating pressure is about 1,3 barg,  the test pressure was set to two times operating pressure or 2.6 barg.

As the results of test were positive at 2.6 barg we decided to push further. We tested the flange up to 8 barg.

Unit passed the test at 8 barg. Achieved result is really good and warrant an optimal sealing and reliable pressure measurements.

Saturday, October 25, 2014

BasicAirData free application available on Google Play™ store

From now on our mobile applications will be available through Google Play™. We've published right now the GPSLogger.

Versions for other devices like Windows Phone and Iphone are not under development; if you are willing to code BasicAirData applications to these devices just enter in contact.

Read our BasicAirData intro to discover on going Open projects, we're looking forward to collaborate with you.

Thursday, October 16, 2014

Pressure test. Front Mounting Flange for 8mm Pitot

A group of users asked for a removable mount for the BasicAirData 8 mm pitot; we decided to implement such a solution. The major design issue was leaking.
Latest design has been just tested up to 8 bar! So finally we can say that the removable airthight and removable flange is ready for be tested on a airframe. As soon as possible we will update the documentation.
Photos from the test

Saturday, October 11, 2014

A state estimation example, GPS aided Altimeter. Part three

Previous post here

In this post we will examine a one degree of freedom inertial unit aided by a GPS altitude measurement. Download from this links the example Matlab or Scilab script.
In the last two articles a basic inertial unit has been presented, in this post we will use measurements from INS and GPS to obtain an altitude sensor. 
Let's pretend GPS sample rate \(T\) is 1 second.  The new altitude measurement equation is
with   \(k\)  positive integer

Between two successive samples the system works exactly as per precedent article equations, the altitude error grows fast.
Every instant \(t=kT\) the fresh measurement arrives from GPS, then a error signal is calculated subtracting the INS position estimate from the GPS measurement. 
$$\delta y_k=\tilde{z}_k-\hat{z}_k=\boldsymbol{C}\delta\hat{x}_k +\nu_z(k) $$
The observability matrix is eye(3,3); The rank is three hence the system is completely observable.
State error evolution is described by the following discrete linear system
$$\boldsymbol{\Phi}=exp(\boldsymbol{A} T)$$
$$\delta \boldsymbol{\hat{x}}^-_{k+1}=\boldsymbol{\Phi}\delta\boldsymbol{\hat{x}}^+_k$$
We set \(var(\nu_z(kT))= 4 m^2\) and \(L=[0.4,0.04,0.002]'\)
Refer to script file for the other parameters values.

The altitude error behavior is plotted in figure 1, the response dynamic is due to the choice of L. 
Figure 1 Error Deviation Plot

With the use of both sensors the altitude deviation is bounded and its final value is near 1.2m; that value is nicely lower than GPS only deviation value that is 2m. To get a stable altitude value are necessary around 50 samples, with our sample rate it is equivalent to 50s. To shape the time response it is necessary, for example, to fine tune the state observer behavior with the poles placement technique.

The filter layout of figure 2 includes an observer and measurements from both the sensors,  In figure 3 you find the digital implementation that includes the INS model; a zero order realization is assumed.

Figure 2 Filter layout
Figure 3 Filter, full digital implementation T is the discrete INS sample time

The shown discretization of INS fits well the system. In figure 4 the response of continuous system is compared with analytical zero order hold discretization and numerical Tustin discretization.
Figure 4 INS only model Input Vs Output and magnified detail view

The figure 2 layout can be used with different kind of sensors;for example,but not only, barometric altimeters, sonar range finders and LIDAR. If the filter is used for altitude measurement the acceleration component along gravity vector and GPS altitude should be used as the system inputs.

The provided zero order hold implementation have a low computational cost and can be implemented on mobile devices. The overall error variance has been bound and its value is lower than that of single instruments.

To the following part

Saturday, September 27, 2014

A state estimation example: GPS aided Altimeter. Part two.

In the previous article  we began the description of a telemetry application for altitude measurement.
We've briefly introduced LTI systems and state observers.

In this post we will examine a one degree of freedom inertial unit. With such a unit it will be possible to directly measure altitude variations. Accompanying scripts for Matlab and Scilab are provided as downloads.
The use of inertial units requires considering the mounting position, as explained in this article.

Let's pretend we want to measure the altitude variation of our aircraft using the measurements from an accelerometer. To obtain an initial value, we define an arbitrary, zero-reference position and an altitude direction. The altitude value from the reference is \(z(t)\) , velocity is \(v(t)\) and acceleration \(a(t)\). From now on, assume that the acceleration measurements are always aligned with \(z\)-axis, which is as if we modelled a simple flat fall but this hypothesis allows us to focus on a specific problem. In a real-word case, other measurements from the rest of the axes of the inertial unit should be used to get the correct acceleration value.

Figure 1 Error variance Vs time for example system, initial variance set to zero
The kinematic equations of our system are:
The measurement system cannot be considered deterministic, since it is prone to measurement errors. We model the accelerometer output as:
\(b(t)\) term represents the random walk of the sensor (Paragraph 3.2.2) and \(\mu_a(t)\) is a Gaussian white noise with given variance \(\sigma^2_{\mu_a}\). The latter term aggregates different kinds of uncertainty sources such as thermal noise.

Random walk noise is defined as \(\dot{b}(t)=\omega_b(t)\) where \(\omega_b(t)\) is Gaussian with given  \(\sigma^2_{\omega_b}\). The average value of \(b(t)\) is zero.

It is interesting to note that while a constant accelerometer error \(\epsilon\) (Paragraph 4.2.1) leads to an error term equal to \(\epsilon t^2/2\), that constant error can be eliminated with an accurate calibration, hence it is not that interesting to take into account and examine.

If we wanted to reconstruct the aircraft altitude based solely on the accelerometer readings and the zero initial altitude assumption, the observed altitude would be:
$$\dot{\hat{z}}(t) = \hat{v}(t)$$
$$\dot{\hat{v}}(t) = \hat{a}(t)$$
Navigation equations are [1][2]
$$\boldsymbol{x}=\begin{bmatrix} \hat{z}\\ \hat{v}\\ \hat{b} \end{bmatrix}$$
As we are interested in error dynamic  we define a new error state variable \(\delta\boldsymbol{x}(t)=\boldsymbol{x}-\boldsymbol{\hat{x}}\)
$$\delta\boldsymbol{x}(t)=\boldsymbol{A}\delta x(t) +\boldsymbol{B}\begin{bmatrix}\mu_a \\ \omega_b \end{bmatrix}$$
General definition of process noise covariance for a continuous system is \(cov(\boldsymbol{\omega}(t), \boldsymbol{\omega}(\tau))=\boldsymbol{Q}(t)\delta(t-\tau)\)
Defining \(\boldsymbol{\Phi}(\tau)=exp(A\tau)\)
The state error covariance is [2]
$$\boldsymbol{P}_x(t)=\boldsymbol{\Phi}(t)\boldsymbol{P}_x(0)\Phi^T(t)+\int_0^{\tau} \boldsymbol{\Phi}(t)\boldsymbol{B}\boldsymbol{Q}\boldsymbol{B}^T\Phi^T(t)d\tau$$
This is one of the rare cases where solution can be calculated in a closed form, terms with subscript 0 denotes initial value.
Let's set \(\sigma_{\mu_a}\)  or velocity random walk to a value of \(5e-3 \frac{m/s}{\sqrt{s}}\) ,  \(\sigma^2_{\mu_a}=2.5e-5\frac{m^2}{s^3}\)  and \(\sigma^2_{\omega_b}=1e-6\frac{m^2}{s^5}\).

As expected and plotted in Figure 1, the use of an inertial sensor alone, leads to an unacceptable altitude error in a few seconds. In the next post this inertial unit will be coupled with a GPS to provide an enhanced measurement. The primary objective is to constrain the error value within bound limits.
1 A.Leon-Garcia (1994), Probability and Random Processes for Electrical Engineering, Prentice Hall
2. Jay A.Farrell (2008), Aided Navigation, Mc Graw Hill

Other readings

Mohinder S. Grewal, Angus P. Andrews (2008), Kalman Fitlering Theory and Practice Using MATLAB, 3rd ed., Wiley, pp 100-101

Friday, September 12, 2014

Terrific Vintage Glider meeting scheduled on september 21th in Cremona, Italy

Terrific Vintage Glider meeting scheduled on september 21th in Cremona, Italy. Basic Air Data will be on the field! 

Thursday, August 7, 2014

New Release, BasicAirData GPSLogger for Android

We have just released a new update of the BasicAirData GPS Logger For Android.
GPS Logger is a simple App for taking track of your position and path, it features:
  • Visualization of GPS position, altitude and speed
  • Recording of tracks in KML and GPX format; you can watch them on Google Earth (also on the smartphone) and on Openstreetmap. Files are saved into /GPSLogger folder on your phone
  • You can use tracks recorded in Openstreetmap editor
  • Compatible with Android 2.2+
  • Battery Frendly (it can record continuously for hours)
In the last article we have already described the features above; this article describes the main changes and the new features included in this update.

Figure 1.The new GPS Logger interface

New functions highlights
  • New and revamped interface, using the ActionBar widgets and the modern Holo Dark Android theme
  • Automatic altitude correction, based on NGA EGM96 Earth Geoid Model
  • Possibility to add Placemarks, with custom label, also while logging tracks
  • More customizable settings for visualization of Altitude and Coordinates   
Let's focus on two of these features. 

Automatic altitude Correction

The main and most interesting feature on which we focused is the correction of GPS altitude error. Recalling the previous article in this link, measurement errors are due to the difference between the earth mean ellipsoid (the zero-reference of the GPS) and the real earth geoid.

There are many models that describe that geoid; the chosen one is the NGA/NASA EGM96, N=M=360 Earth Gravitational Model, described by the binary geoid height file hosted on National Geospatial-Intelligence Agency site.

There are newer and more detailed geoid models than EGM96, but this one is accurate enough for smartphones; it contains a reasonable amount of freely loadable data, so this is the best suitable for Smartphone usage.

When The Binary EGM96 Automatic Correction is selected for the first time into the Setting Screen(See Figure 2) the geoid Height file is downloaded from the NGA/NASA website. The size of that file is about 2 MB.

Figure 2.  Enabling EGM96 Correction for the first time

The software code is tailored for the format of this file. The file you downloaded is an unformatted direct access file. The data is arranged in 721 arrays of 1440 records each. Arrays are ordered from north to south. Records are disposed from west to east starting at the prime meridian (0 E) and ending 15 arc-minutes west of the prime meridian (359.75 E). Data unit of measurement is the centimeter. The grid of 721 x 1440 values represents the distance from Ellipsoid and Earth Geoid, or geoid deviation, every 15 arc-minutes.

At every position fix, by means of a bi-linear interpolation, the application recalculates the correction value for the actual position.

To have a general idea of the correction dynamic we did some basic analysis of geoid heights file. The graph in figure 3 shows the output of this Scilab script (a parser that calculates the maximum gradient value of the whole grid) WW15MGH_max_grad.sce, that represents the area with the maximum gradient of the EGM96 grid. That place (approx. 30 N, 80 E, situated on Himalaya Mountains) has a gradient value of:
$$\mid\nabla\boldsymbol{h}\mid=974 cm / 15 arc-second = 9.74 m / 24000 m = 4.058e-4$$
That means than correction value may change at least 1 cm every 24 m.

Figure 3. Area of maximum gradient of EGM96 grid

Let's pretend that we're logging with GPS Logger at 1 Hz on a car running at 90 km/h. Considering the previous gradient data then we should expect that our EGM96 correction value might change at every single GPS sample.

The application minimum correction value is then chose to be constant and equal to 1 cm, even if a coarse discretized GPS measurement is used. In such a way the overall measurement uncertainty is minimized in a wide set of operating conditions.

Placemark creation

The last feature we describe here is the placemark creation.
It is accessible from the placemark icon on the actionbar, it allows to mark the succeeding fix with a mark symbol upon map and to specify a label that describes it. That function is useful to remember special locations during your path, or to write some information during mapping (by instance for Openstreetmap contribution). Alternatively you can dictate the description with Google Vocal Dictation.
It is possible to record a placemark during track recording, and also in while in pause.

Figure 4 Placemark creation, to add info on OpenStreetMap

 Installation of the app is straightforward, you may find detailed information here.

At a glance, the installation consist in three steps.
  1. Download the package here
  2. Copy the GPSLogger.apk file on your cellphone SD card and
  3. Install it.
Alternatively you can use the QR-CODE here below.

Figure 5. QR-Code for GPS Logger application download

The app have been developed using Android Studio under Fedora 19.
It's 100% free and Open Source. So you can freely download and install it, share it with your friends, and also download the whole package containing the source code in this link. The code is written in Java and commented by the author. If any doubt arise just ask.

Figure 6. A long mountain route. Data from GPS Logger shown in Google Earth

The application is basic but useful in many DIY surveying tasks, it is be proven to be ideal to support some basic calibration procedures for air data instruments.

Thursday, July 31, 2014

A state estimation example: GPS-aided altimeter

The value of telemetry data can be enhanced by the application of basic estimation and modeling techniques. In most real-life cases, even in the simplest systems, it is either too expensive or actually impossible to measure every relevant variable. Thus, one needs to make do with a reduced number of sensors, which are also quite noisy.

In the next few articles, I will present a complete application example in the form of a miniseries. The core idea is to investigate the performance of an algorithm that uses data from accelerometers and a GPS to enhance altitude measurement. That method should be implementable on mobile devices and can enhance the previously implemented altimeter. The very same algorithm can be used to post-process basic flight telemetry data.

Miniseries summary:
  • Introduction post: the state estimation problem and state observer
  • Altitude measurement enhancement with accelerometer, performance evaluation
  • Altitude measurement enhancement with accelerometer and GPS data, performance evaluation

Working code will be posted where needed.
In order to use telemetry data in a meaningful way, it is necessary to have a mathematical model of our physical system. One possible representation is the state space representation.
A continuous linear system is represented as 
$$\boldsymbol{y}(t)= \boldsymbol{C}(t)\boldsymbol{x}(t)+ \boldsymbol{D}(t)\boldsymbol{u}(t)$$
Bold notation indicates vectors or matrices. All quantities are time dependent. This is a fully deterministic model, whose the mathematical representation is a set of ordinary differential equations.
\(\boldsymbol{x}(t)\) is the vector containing the physical system state variables, \(\boldsymbol{y}(t)\) represents the output variables and \(\boldsymbol{u}(t)\) is the system input. All the other matrices describe the relationships between those variables.
However, in our example case, we can make some assumptions towards the simplification of the system:
$$\boldsymbol{\dot{x}}(t)= \boldsymbol{A}\boldsymbol{x}(t)+ \boldsymbol{B}\boldsymbol{u}(t)$$
$$\boldsymbol{y}(t)= \boldsymbol{C}\boldsymbol{x}(t)+ \boldsymbol{D}\boldsymbol{u}(t)$$
Now system matrices are time independent. Last set of equations is a model for a Linear Time Invariant system. If \(\boldsymbol{x}\) has a dimension of \(n \times 1\), then we say that the system has a dimension of \(n\).
To better explain the state space representation I'll use a simple mechanical system. Suppose we need the mathematical model of a falling apple.

Figure 1 Falling Apple

To simplify the math, suppose the falling apple is not rotating. Assume also that the apple is moving through the void, thus aerodynamics are not involved. The apple will fall vertically.

Figure 2 Model definition

The system dynamics can be described using Newton's second law:
$$\boldsymbol{F}=m\boldsymbol{a_z}= m\boldsymbol{g} \rightarrow a_z=g$$
Using standard notation we can write:
$$\begin{bmatrix}\dot{p} \\ \dot{v} \end{bmatrix} =\begin{bmatrix}0 &1 \\ 0& 0 \end{bmatrix}    \begin{bmatrix}p \\ v \end{bmatrix}+ \begin{bmatrix}0 \\ 1 \end{bmatrix} g  $$
Note that we have chosen not to model the output \(\boldsymbol{y}(t)\). It is not necessary in this example.
You can observe that such a model is not completely realistic. According to it, the apple will fall to infinity with a uniformly accelerated motion. As in most occasions the model is accurate only around specific conditions: the previous apple model is accurate for a short fall in the void, before the apple hits the plate. Commonly the model of a phenomenon is not linear but can be considered linear around predefined state variables values. That's one reason why linear models are useful for modeling real-world systems.

Figure 3 Mass-Spring-Damper System

Let’s move over to a more complex system and examine a complete numerical example. The system is composed by a mass, a spring and a damper. It is a classic literature example, for which you can find a description in this link.
Download in these links the Scilab and Matlab script files for the Mass-Spring-Damper example.

Let's suppose that we want to know the position and the speed of the mass with a single position sensor. I know it's simple to find the solution, however the systematic solution that I will show here can be used on more complex systems. in particular it will be used on the inertial navigation system which will be presented later.

State estimation consists of getting the values of state variables without measuring them. We will achieve such a result with the use of a state observer.

The main working hypothesis is that the physical system is correctly modeled, in this case as an LTI system  of known \(\boldsymbol{A}, \boldsymbol{B}, \boldsymbol{C}, \boldsymbol{D}\) matrices. The model of the system is supposed to describe exactly the physical system behavior. We do not need to identify the parameters or the order our system but we are only looking for the state space variables values. That the main idea is to get measurement values from a limited set of sensors and obtain the observed state estimation. To that goal, we will use a Luenberger Observer architecture, an estimator which is a closed loop state observer.

Figure 4 Estimation layout

Estimated variables will be indicated with the circumflex or hat notation.
Returning to the mass-spring-damper system the continuous time estimator will be described by the following system of equations. The only available sensor measurement  is the mass position. The estimator gain is contained in the \(\boldsymbol{L}\) matrix.
$$\boldsymbol{\hat{y}}(t)=\boldsymbol{C}\boldsymbol{\hat{x}}(t)+ \boldsymbol{D}\boldsymbol{u}(t)$$
Given spring constant \(k\), damping constant \(c\) and body mass \(m\), we get
$$\boldsymbol{A}=\begin{bmatrix}0&1 \\ -k/m&–c/m \end{bmatrix}$$
$$\boldsymbol{B}=\begin{bmatrix} 0  \\ 1/m \end{bmatrix}$$
$$\boldsymbol{C}=\begin{bmatrix} 1 & 0\end{bmatrix}$$
$$\boldsymbol{D}=\begin{bmatrix} 0\end{bmatrix}$$
For presentation purposes, in the provided script files, the simulation state starts from [1 ; 0] and the system input is set to zero. You can see how the position oscillates with time and how well the observers track the velocity state.

Figure 5 Time response of system and state observer

Prior to the use of an observer it's necessary to check the observability matrix rank. If the rank is equal to the system order then it is possible to observe every state of the system. For our system \(O=[C;C*A]\) and the rank of \(O\) is two, equal to the rank of our system \(n\) which is also is two. Hence the system is completely observable (command rank(O)). Another property that should be checked is the system stability. For now we will avoid to deal with unstable systems. By issuing the command spec(A) in Scilab or eig(A) in Matlab we get the system eigenvalues, which are also called poles of the system. If the real part of every pole is negative then the system is stable. In our case the eigenvalues are a complex conjugate pair -0.1 ±0.53i, and the system is stable.
With complex eigenvalues the state variables trajectories will oscillate when moving from one stable solution to another. This behavior is evident in the plots of figure 5.

If we define the tracking error as \(\boldsymbol{e}(t)=\boldsymbol{\hat{x}}(t)-\boldsymbol{x}(t)\) we get 
Hence, the error dynamic response is directly defined by \(\boldsymbol{L}\) matrix. The poles of \(\boldsymbol{A-LC}\) should be stable to have a usable estimator. In our example we get the stable complex poles -0.6 ±2.1i, which generate an oscillating, but stable, trajectory that is evident in figure 5.

With the dynamic of the error so well defined, it is possible to shape the response, for example with a pole placement method.  If we want two real poles at -5, the we can issue the Scilab command Lt=(ppol(A',C',[-0.5 ;-0.5]))'. Afterwards we can check the closed loop poles values with spec(A-Lt*C). In Matlab the commands are Lt= place(A',C',p).' and eig(A-Lt*C).

After this little introduction we have the tools to understand the altitude measurement system enhanced with accelerometers and GPS. A key point to our further analysis will be the characterization of the error variance of the state estimation.