Inversion of complex matrices: usefulness and how to obtain in MSL

[ceraolo]

Inversion of complex matrices might seem a too esoteric task to be inserted in MSL.

In my opinion it is not. One part of MSL, Electrical.Analog.QuasiStationary relies completely on complex numbers. I've recently found a very important application of Quasi Stationary library: phasor simulation of an AC network which interferes with the dynamic systems constituted by trains: the feeding system of High-Speed railway lines. To do this satisfactorily, complex inversion is needed. I omit here details for brevity, but in case complex inversion is added to MSL I can supply a significant, yet sufficiently simple, as an exampleto be possibly included in MSL. For this application the best solution would be a "cinv" Modelica function which directly inverts a complex matrix. In fact: inversion is only needed at a system initial set-up, to create the parameter complex matrix of this system so

• there is no need of high speed computation
• inversion is to be made in a function, i.e. in an algorithm
• the matrix to be inverted is expected to be always well-conditioned.

Naturally I have several options to solve my specific problem , e.g.:

• imposing the product of the matrix to be inverted times the expected inversion to be equal to identity in an equation section of a model
• making my own C code for inversion

But I prefer to work on a more general solution which could imply an enhancement of MSL's usability.

It seems to me that adding complex inverse should not be a strong task, given the existing Lapack functions (cgetrf() , cgetri()). Since currently Modelica interfaces to C programs through real numbers, a simple Modelica wrapper would be needed to extract the two real matrix from the complex matrix to be inverted, send them to a C wrapper which combines them back into a complex matrix, use Lapack for inversion and make the inverse path to give the Modelica caller function the result.

I could to this myself and propose the result for possible insertion in MSL. I only need to understand how to link cgetri() and the corresponding lapack dependencies to my Modelica installation. I repeat, it f there is interest, I can share everything I do on this ticket's topic to the Modelica community.

Can someone help?

[beutlich]

Yes, I can help for the technical stuff, but I cannot answer the questions on usefulness and interest.

[ceraolo]

Yes, I can help for the technical stuff, but I cannot answer the questions on usefulness and interest. Naturally I think this question is very interesting and useful, but this is just my opinion. I must find a sound way to convince others. I plan at least one paper showing advantages of usage of phasors (i.e. complex numbers) for railway simulations, and can try to do this using workarounds If I succeed, I can share the paper to see if this raises interest. I think QuasiStationary library of MSL has not been used much up to know. An indirect demonstration is that my attempts have generated some tickets and comments on OpenModelica track, on this forum, and even a couple of reports to Modelon regarding Dymola.
This is because time domain is much more natural. But in systems in which me must follow trains for minutes or fractions of hours, following 50 Hz sines is tremendously demanding of computer power, and largely unnecessary. And since the railway supply systems are single phase, this cannot be solved with Park's transformation, as is brilliantly done in Power System library (there this is possible because they are three-phase systems).

So, to conclude, I propose to keep this ticket a bit on hold to see whether I'me able to go forward with workarounds. In case I have success, maybe I could convince many of the usefulness of empowering Modelica's capability of working with complex numbers; in case I'm really stuck I could ask here again at least for a limited help to go forward.

[christiankral]

The inverse of a (complex) matrix is something like a basic functionality. Therefore, on my opinion it makes very much sense to add such a function to the MSL, as the inverse function also exists in Modelica.Math.Matrices.inv for real variables. The function shall be placed in Modelica.ComplexMath.Matrices.inv to keep the structure symmetric with the real function.

I am in favor of implementing such inverse function.

[AHaumer]

I am also in favor. Maybe somebody has the knowledge or time to search an appropriate algorithm.

[sjoelund]

http://www.netlib.org/lapack/lapack-3.1.1/html/zgetri.f.html - but the Modelica external interface does not include complex numbers and calling Fortran functions from external "C" is rather complicated.

[ceraolo]

In my original post I mentioned C Lapack routines, and a technique to overcome the limitation that Modelica external interface does not include complex numbers

[sjoelund]

@ceraolo the problem is that it seems most Modelica tools use different approaches to the lapack and blas libraries. It's not as simple as calling a routine from an external C library. Yes, it's the only thing needed for your proposal to work, but you would possibly need all tool vendors to agree on using for example clapack (rather than Fortran like MSL uses) as its interface (also, some tools might have included only the relevant double-precision functions for simplicity). It could be a good change to use clapack and only have external "C" in MSL and possibly deprecate Fortran from the specification; it would make implementation of a Modelica compiler simpler...

[ceraolo]

aha! I suspected that I missed something. And make all (or nearly all) tool vendors to agree on a common library choice is rather complicated. For what I can see from my standpoint the fact that "most Modelica tools use different approaches to the lapack and blas libraries" is a weakness of Modelica. As you say, in case it is agreed to switch from Fortran to C libraries for built-in modelica functions, could give advantages that go well beyond this ticket's scope. But the advantages must be compared with the effort needed to do the shift, obviously. Ehm... I thought that clapack kept pace with standard (Fortran) Lapack. Looking at the Internet is seems that, instead, that it is not updated since 10 years ago. It this is true, this is an additional problem for the idea of switching everything to "C" in MSL. Maybe using the C language APIs of standard (Fortran) Lapack is ok?

[max-privato]

In the last row of my last message I referred to "Standard C language APIs for LAPACK", as in http://www.netlib.org/lapack/#_presentation and http://www.netlib.org/lapack/lapacke.html

[HansOlsson]

Apart from the problem with link with Lapack etc there are two additional comments related to this:

Inverting a matrix is often not a good idea, and solving linear systems of equations is often better. Obviously we can still make inversion available - but only having inversion without factoring and solving linear systems is not a good idea. Thus it is not a matter of just adding one function, but several in a coherent way.

There are also two ways of handling complex matrices - either as matrix of complex numbers or as two real matrices. Lapack uses the first and that seems natural in Modelica as well, and it seems Modelica's C-interface should handle it as well. (But not for Fortran.)

[max-privato]

reply to HansOlsson

Inverting a matrix is often not a good idea, and solving linear systems of equations is often better.

I agree. there are however cases in which inversion is the best solution. In my ticket's description I mentioned that I could supply details, and I try to do this now. Consider a case in which pre-processing of user's data requires the inversion of complex matrices. This can be done in an initial equation section of even using binding equations calling cinv(). In case the matrix to be inverted is well-conditioned, this is good and fast: only one inversion before the simulation starts is needed.

Obviously one can find the inverse MI of M asking a modelica tool to solve MI x M=I where I is the identity matrix, and this is what I currently do. However, having the possibility of directly inverting the matrix would make software cleaner and easier to use, in cases (as mine) in which it requires several inversions while preparing data for simulations.

Naturally I'm not asking a programming effort just for a problem of mine. I suppose that the usage of matrix inversion when pre-processing (i.e. when determining parameters needed for simulations which are a complicated function of the input parameters) is a rather general situation. If I understand well, Christian Kral agrees on this (see his comment). After all, the inversion of a real matrix is already available, and the objection to the usage of matrix inversion is the same for real and complex ones.

Obviously we can still make inversion available - but only having inversion without factoring and solving linear systems is not a good idea. Thus it is not a matter of just adding one function, but several in a coherent way.

I'm not sure. If you mean that to fix this ticket a full set of complex functions is to be created, I wonder whether this requires too much effort. I really don't know. Maybe just inversion is just what most people might need?

@christiankral can you comment on this?

[tobolar]

There is a possible candidate for implementation realized in the Modelica_LinearSystems2 library: Modelica_LinearSystems2.WorkInProgress.Math.Complex.Matrices.inv

[HansOlsson]

I'm not sure. If you mean that to fix this ticket a full set of complex functions is to be created, I wonder whether this requires too much effort. I really don't know. Maybe just inversion is just what most people might need?

I doubt that - and even if it is the case for "most" there will be others. They will likely use what is available if it solves the problem, even if not ideal - and then they will complain that the inverse generates bad results.

Note that I'm not saying that we shouldn't add inverse for complex matrices - just that we also add solve-functions.

[max-privato]

There is a possible candidate for implementation realized in the Modelica_LinearSystems2 library: Modelica_LinearSystems2.WorkInProgress.Math.Complex.Matrices.inv

Sounds good. I'll look at this soon.

Wisely enough, the function comment says: "Inverse of a comlex matrix (try to avoid, use function solve(..) instead)" So users are warned that inv() usage must be taken with care. Still, there are cases in which this function results useful.

[christiankral]

The solution of complex equations is usually related with of initial conditions. So we have some options, as I understand the discussion:

1. Solve of a linear system of complex equations by writing complex Modelica equations in matrix form; this is what we can already do today; this is basically the way, how @max-privato solves his issue currently; similar to option 2. we may run into accuracy problems in this case, if the linear system of complex equations gets too big or bad conditioned: @HansOlsson do I understand this correct?
2. Implement a complex inverse function in Modelica; as I understand the discussion such a function may cause bad results depending on the size and conditioning of such a complex value problem
3. Implement a solve function which includes factoring to get reliable results of solving a bad conditioned or numerically demanding linear system of complex values;

Until today we use option 1. since options 2. and 3. are not yet available in the MSL. I think it makes sense to include option 2. and 3. to the MSL. In some simple cases, option 2 may be sufficient for simplicity reasons. In case of more complex and numerically demanding systems, option 3. may even be required to solve bad conditioned systems of equations. So this may not be a luxury problem but rather required for some kinds of problems.

[christiankral]

@tobolar As of Modelica_LinearSystems2 I do not understand if it is planned to include this library to the MSL or not. Having a solution in the MSL looks quite beneficial to me.

[HansOlsson]

@christiankral Regarding 1: If we write a complex matrix equation and have tools solve this (using linear algebra with real matrices) that should not cause any accuracy problems. However, as I understand it might be slightly less convenient to write it in this form.

Regarding the other two options: If we use a complex inverse function it will usually internally factor the matrix first and then solve (or use a slightly improve solve-function). Thus exposing a solve function doesn't bring in more code - it just allows us to call it more efficiently.

[max-privato]

@tobolar I tried a simple example with Dymola (which is declared to be the tool under which the library is tested) but I got: dymosim started ... "dsin.txt" loading (dymosim input file) (j>=1)&&(j<=dim) The following error was detected at time: 0 Index out of bounds The stack of functions is: Modelica_LinearSystems2.WorkInProgress.Math.Complex.Matrices.inv( ZrInv) Error: Failed to start model.

This is strange since the only parameter I have to supply is the matrix to be inverted. BTW It was not clear to me where to locate zlapack.ib. I used the one in https://svn.modelica.org/projects/Modelica_EmbeddedSystems/branches/Modelica_LinearSystems2-with-more-integrators/Modelica_LinearSystems2/

Do you know if there where is to be found the correct zlapack.lib and if it is a 32 or 64 bit library? In addition, do you know is there is somewhere a working inv() example?

[tobolar]

@christiankral

I do not understand if it is planned to include this library to the MSL or not.

I guess there is a wish to include the libraray into the MSL, but no plan on how to proceed. So moving reasonable functionalities into the MSL step-by-step could be a feasible way.