Solving a system of equations using Python/Scipy for a set of measurements

Question

I have an physical instrument of measurement (force platform with load cells) which gives me three values, A, B and C. It happens, though, that these values - that should be orthogonal - actually are somewhat coupled, due to physical characteristics of the measuring device, which causes cross-talk between applied and returned values of force and torque.

Then, it is recommended that a calibration matrix be used to transform the measured values into a better estimate of the actual values, like this:

The problem is that it is necessary to perform a SET of measurements, so that different measured(Fz, Mx, My) and actual(Fz, Mx, My) are least-squared to get some C matrix that works best for the system as a whole.

I can solve Ax = B problems with scipy.linalg.lststq, or even scipy.linalg.solve (giving an exact solution) for ONE measurement, but how should I proceed to consider a set of different measurements, each one with its own equation giving a potentially different 3x3 matrix?

Any help is much appreciated, thanks for reading.

Answer 1

I posted a similar question containing just the mathematical part of this at math.stackexchange.com, and this answer solved the problem:

math.stackexchange.com/a/232124/27435

In case anyone have a similar problem in the future, here is the almost literal Scipy implementation of that answer (first lines are initialization boilerplate code):

import numpy
import scipy.linalg

### Origin of the coordinate system: upper left corner!
"""
    1----------2
    |          |
    |          |
    4----------3
"""

platform_width = 600
platform_height = 400

# positions of each load cell (one per corner)
loadcell_positions = numpy.array([[0, 0],
                                  [platform_width, 0],
                                  [platform_width, platform_height],
                                  [0, platform_height]])

platform_origin = numpy.array([platform_width, platform_height]) * 0.5

# applying a known force at known positions and taking the measurements
measurements_per_axis = 5
total_load = 50

results = []
for x in numpy.linspace(0, platform_width, measurements_per_axis):
    for y in numpy.linspace(0, platform_height, measurements_per_axis):
        position = numpy.array([x,y])
        for loadpos in loadcell_positions:
            moments = platform_origin-loadpos * total_load
            load = numpy.array([total_load])
            result = numpy.hstack([load, moments])
            results.append(result)
results = numpy.array(results)
noise = numpy.random.rand(*results.shape) - 0.5
measurements = results + noise

# now expand ("stuff") the 3x3 matrix to get a linearly independent 3x3 matrix
expands = []
for n in xrange(measurements.shape[0]):
    k = results[n,:]
    m = measurements[n,:]

    expand = numpy.zeros((3,9))
    expand[0,0:3] = m
    expand[1,3:6] = m
    expand[2,6:9] = m
    expands.append(expand)
expands = numpy.vstack(expands)

# perform the actual regression
C = scipy.linalg.lstsq(expands, measurements.reshape((-1,1)))
C = numpy.array(C[0]).reshape((3,3))

# the result with pure noise (not actual coupling) should be
# very close to a 3x3 identity matrix (and is!)
print C

Hope this helps someone!