Proposition to create a complex matrix class

iidm/iidm-extensions/src/main/java/com/powsybl/iidm/network/extensions/util/FortescueUtil.java

@@ -7,8 +7,11 @@
  */
 package com.powsybl.iidm.network.extensions.util;
 
+import com.powsybl.math.matrix.ComplexMatrix;
 import com.powsybl.math.matrix.DenseMatrix;
+import org.apache.commons.math3.complex.Complex;
 import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
+import org.apache.commons.math3.util.FastMath;
 import org.apache.commons.math3.util.Pair;
 
 /**
@@ -17,76 +20,69 @@
 public final class FortescueUtil {
 
     public enum SequenceType {
-        POSITIVE,
-        NEGATIVE,
-        ZERO
+        POSITIVE(1),
+        NEGATIVE(2),
+        ZERO(0);
+
+        private final int num;
+
+        SequenceType(int num) {
+            this.num = num;
+        }
+
+        public int getNum() {
+            return num;
+        }
     }
 
     public static DenseMatrix getFortescueMatrix() {
 
         // [GA]   [ 1  1  1 ]   [G0]
         // [GB] = [ 1  a²  a] * [G1]
         // [GC]   [ 1  a  a²]   [G2]
-        return getFortescueOrInverseMatrix(false);
+        return createComplexMatrix(false).toRealCartesianMatrix();
     }
 
     public static DenseMatrix getFortescueInverseMatrix() {
 
         // [G0]       [ 1  1  1 ]   [GA]
         // [G1] = 1/3 [ 1  a  a²] * [GB]
         // [G2]       [ 1  a² a ]   [GC]
-        return getFortescueOrInverseMatrix(true);
+        return createComplexMatrix(true).toRealCartesianMatrix();
     }
 
-    public static DenseMatrix getFortescueOrInverseMatrix(boolean isInverse) {
-        DenseMatrix mFortescueOrInverse = new DenseMatrix(6, 6);
+    public static ComplexMatrix createComplexMatrix(boolean inverse) {
+        // [G1]   [ 1  1  1 ]   [Gh]
+        // [G2] = [ 1  a²  a] * [Gd]
+        // [G3]   [ 1  a  a²]   [Gi]
+
+        Complex a = new Complex(-0.5, FastMath.sqrt(3.) / 2);
+        Complex a2 = a.multiply(a);
 
-        double t = 1;
-        double signA = 1;
-        if (isInverse) {
+        double t = 1.;
+        Complex c1 = a;
+        Complex c2 = a2;
+        if (inverse) {
             t = 1. / 3.;
-            signA = -1;
+            c1 = a2.multiply(t);
+            c2 = a.multiply(t);
         }
+        Complex unit = new Complex(t, 0);
 
-        //column 1
-        mFortescueOrInverse.add(0, 0, t);
-        mFortescueOrInverse.add(1, 1, t);
-
-        mFortescueOrInverse.add(2, 0, t);
-        mFortescueOrInverse.add(3, 1, t);
-
-        mFortescueOrInverse.add(4, 0, t);
-        mFortescueOrInverse.add(5, 1, t);
-
-        //column 2
-        mFortescueOrInverse.add(0, 2, t);
-        mFortescueOrInverse.add(1, 3, t);
-
-        mFortescueOrInverse.add(2, 2, -t / 2.);
-        mFortescueOrInverse.add(2, 3, signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(3, 2, -signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(3, 3, -t / 2.);
-
-        mFortescueOrInverse.add(4, 2, -t / 2.);
-        mFortescueOrInverse.add(4, 3, -signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(5, 2, signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(5, 3, -t / 2.);
-
-        //column 3
-        mFortescueOrInverse.add(0, 4, t);
-        mFortescueOrInverse.add(1, 5, t);
+        ComplexMatrix complexMatrix = new ComplexMatrix(3, 3);
+        complexMatrix.set(0, 0, unit);
+        complexMatrix.set(0, 1, unit);
+        complexMatrix.set(0, 2, unit);
 
-        mFortescueOrInverse.add(2, 4, -t / 2.);
-        mFortescueOrInverse.add(2, 5, -signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(3, 4, signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(3, 5, -t / 2.);
+        complexMatrix.set(1, 0, unit);
+        complexMatrix.set(1, 1, c2);
+        complexMatrix.set(1, 2, c1);
 
-        mFortescueOrInverse.add(4, 4, -t / 2.);
-        mFortescueOrInverse.add(4, 5, signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(5, 4, -signA * t * Math.sqrt(3.) / 2.);
-        mFortescueOrInverse.add(5, 5, -t / 2.);
+        complexMatrix.set(2, 0, unit);
+        complexMatrix.set(2, 1, c1);
+        complexMatrix.set(2, 2, c2);
 
-        return mFortescueOrInverse;
+        return complexMatrix;
     }
 
     public static Vector2D getCartesianFromPolar(double magnitude, double angle) {

math/pom.xml

@@ -40,6 +40,10 @@
 
     <dependencies>
         <!-- Compilation dependencies -->
+        <dependency>
+            <groupId>org.apache.commons</groupId>
+            <artifactId>commons-math3</artifactId>
+        </dependency>
         <dependency>
             <groupId>com.google.guava</groupId>
             <artifactId>guava</artifactId>

math/src/main/java/com/powsybl/math/matrix/ComplexMatrix.java

@@ -0,0 +1,129 @@
+/**
+ * Copyright (c) 2022, RTE (http://www.rte-france.com)
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+package com.powsybl.math.matrix;
+
+import org.apache.commons.math3.complex.Complex;
+import org.apache.commons.math3.util.FastMath;
+
+import java.util.Objects;
+
+/**
+ * @author Jean-Baptiste Heyberger {@literal <jbheyberger at gmail.com>}
+ */
+public class ComplexMatrix {
+
+    private static final double EPSILON = 0.00000001;
+
+    private final DenseMatrix realPartMatrix;
+    private final DenseMatrix imagPartMatrix;
+
+    public ComplexMatrix(int rowCount, int columnCount) {
+        realPartMatrix = new DenseMatrix(rowCount, columnCount);
+        imagPartMatrix = new DenseMatrix(rowCount, columnCount);
+    }
+
+    public void set(int i, int j, Complex complex) {
+        Objects.requireNonNull(complex);
+        realPartMatrix.set(i, j, complex.getReal());
+        imagPartMatrix.set(i, j, complex.getImaginary());
+    }
+
+    public int getRowCount() {
+        return realPartMatrix.getRowCount();
+    }
+
+    public int getColumnCount() {
+        return realPartMatrix.getColumnCount();
+    }
+
+    public Complex get(int i, int j) {
+        return new Complex(realPartMatrix.get(i, j), imagPartMatrix.get(i, j));
+    }
+
+    public static ComplexMatrix createIdentity(int nbRow) {
+        ComplexMatrix complexMatrix = new ComplexMatrix(nbRow, nbRow);
+        for (int i = 0; i < nbRow; i++) {
+            complexMatrix.realPartMatrix.set(i, i, 1.);
+        }
+        return complexMatrix;
+    }
+
+    public ComplexMatrix transpose() {
+        ComplexMatrix transposed = new ComplexMatrix(getColumnCount(), getRowCount());
+        for (int i = 0; i < getColumnCount(); i++) {
+            for (int j = 0; j < getRowCount(); j++) {
+                transposed.realPartMatrix.set(i, j, realPartMatrix.get(j, i));
+                transposed.imagPartMatrix.set(i, j, imagPartMatrix.get(j, i));
+            }
+        }
+        return transposed;
+    }
+
+    public ComplexMatrix scale(Complex factor) {
+        Objects.requireNonNull(factor);
+        ComplexMatrix scaled = new ComplexMatrix(getRowCount(), getColumnCount());
+        for (int i = 0; i < getRowCount(); i++) {
+            for (int j = 0; j < getColumnCount(); j++) {
+                scaled.realPartMatrix.set(i, j, realPartMatrix.get(i, j) * factor.getReal() - imagPartMatrix.get(i, j) * factor.getImaginary());
+                scaled.imagPartMatrix.set(i, j, imagPartMatrix.get(i, j) * factor.getReal() + realPartMatrix.get(i, j) * factor.getImaginary());
+            }
+        }
+        return scaled;
+    }
+
+    public ComplexMatrix scale(double factor) {
+        return scale(new Complex(factor, 0.));
+    }
+
+    // utils to switch between complex and real cartesian representation of a complex matrix
+    public DenseMatrix toRealCartesianMatrix() {
+        DenseMatrix realMatrix = new DenseMatrix(getRowCount() * 2, getColumnCount() * 2);
+        for (int i = 0; i < getRowCount(); i++) {
+            for (int j = 0; j < getColumnCount(); j++) {
+                Complex complexTerm = new Complex(realPartMatrix.get(i, j), imagPartMatrix.get(i, j));
+                realMatrix.add(2 * i, 2 * j, complexTerm.getReal());
+                realMatrix.add(2 * i + 1, 2 * j + 1, complexTerm.getReal());
+                realMatrix.add(2 * i, 2 * j + 1, -complexTerm.getImaginary());
+                realMatrix.add(2 * i + 1, 2 * j, complexTerm.getImaginary());
+            }
+        }
+
+        return realMatrix;
+    }
+
+    public static ComplexMatrix fromRealCartesian(DenseMatrix realMatrix) {
+        Objects.requireNonNull(realMatrix);
+        int columnCount = realMatrix.getColumnCount();
+        int rowCount = realMatrix.getRowCount();
+        if (columnCount % 2 != 0 || rowCount % 2 != 0) { // dimensions have to be even
+            throw new MatrixException("Incompatible matrices dimensions to build a complex matrix from a real cartesian");
+        }
+
+        ComplexMatrix complexMatrix = new ComplexMatrix(rowCount / 2, columnCount / 2);
+        for (int i = 0; i < rowCount / 2; i++) {
+            for (int j = 0; j < columnCount / 2; j++) {
+
+                int rowIndexInCartesian = 2 * i;
+                int colIndexInCartesian = 2 * j;
+
+                // Before building the complex matrix term, check that the 4x4 cartesian bloc can be transformed into a Complex term
+                double t11 = realMatrix.get(rowIndexInCartesian, colIndexInCartesian);
+                double t12 = realMatrix.get(rowIndexInCartesian, colIndexInCartesian + 1);
+                double t21 = realMatrix.get(rowIndexInCartesian + 1, colIndexInCartesian);
+                double t22 = realMatrix.get(rowIndexInCartesian + 1, colIndexInCartesian + 1);
+
+                if (FastMath.abs(t11 - t22) > EPSILON || FastMath.abs(t12 + t21) > EPSILON) {
+                    throw new MatrixException("Incompatible bloc matrices terms to build a complex matrix from a real cartesian");
+                }
+
+                complexMatrix.set(i, j, new Complex(t11, t21));
+            }
+        }
+
+        return complexMatrix;
+    }
+}

math/src/test/java/com/powsybl/math/matrix/ComplexMatrixTest.java

@@ -0,0 +1,56 @@
+/**
+ * Copyright (c) 2022, RTE (http://www.rte-france.com)
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+package com.powsybl.math.matrix;
+
+import org.apache.commons.math3.complex.Complex;
+import org.junit.jupiter.api.Test;
+
+import static org.junit.jupiter.api.Assertions.assertEquals;
+
+/**
+ * @author Jean-Baptiste Heyberger {@literal <jbheyberger at gmail.com>}
+ */
+class ComplexMatrixTest {
+
+    @Test
+    void complexMatrixTest() {
+
+        ComplexMatrix cm = new ComplexMatrix(2, 2);
+        cm.set(0, 0, new Complex(1., 2.));
+        cm.set(0, 1, new Complex(3., 4.));
+        cm.set(1, 0, new Complex(5., 6.));
+        cm.set(1, 1, new Complex(7., 8.));
+
+        // create the real equivalent matrix of complex matrix
+        DenseMatrix rm = cm.toRealCartesianMatrix();
+
+        assertEquals(cm.get(0, 0).getReal(), rm.get(0, 0));
+        assertEquals(cm.get(1, 1).getImaginary(), -rm.get(2, 3));
+
+        ComplexMatrix cmFromReal = ComplexMatrix.fromRealCartesian(rm);
+
+        for (int i = 0; i < cm.getRowCount(); i++) {
+            for (int j = 0; j < cm.getColumnCount(); j++) {
+                assertEquals(cm.get(i, j), cmFromReal.get(i, j));
+            }
+        }
+
+        ComplexMatrix tm = cm.transpose();
+        assertEquals(cm.get(0, 1), tm.get(1, 0));
+
+        ComplexMatrix im = ComplexMatrix.createIdentity(2);
+        assertEquals(im.get(0, 1), new Complex(0, 0));
+        assertEquals(im.get(1, 1), new Complex(1, 0));
+
+        ComplexMatrix sm = cm.scale(new Complex(2., 1.)); //ComplexMatrix.getMatrixScaled(cm, new Complex(2., 1.));
+        assertEquals(sm.get(0, 0), new Complex(0., 5.));
+
+        sm = cm.scale(3.);
+        assertEquals(sm.get(0, 0), new Complex(3., 6.));
+
+    }
+}