Skip to content

Commit

Permalink
Proposition to create a complex matrix class (#2835)
Browse files Browse the repository at this point in the history
Signed-off-by: JB-H <jbheyberger@gmail.com>
Co-authored-by: Geoffroy Jamgotchian <geoffroy.jamgotchian@rte-france.com>
  • Loading branch information
JB-H and geofjamg authored May 7, 2024
1 parent 2f17f08 commit f8ef9c6
Show file tree
Hide file tree
Showing 4 changed files with 232 additions and 47 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -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;

/**
Expand All @@ -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) {
Expand Down
4 changes: 4 additions & 0 deletions math/pom.xml
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
129 changes: 129 additions & 0 deletions math/src/main/java/com/powsybl/math/matrix/ComplexMatrix.java
Original file line number Diff line number Diff line change
@@ -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;
}
}
56 changes: 56 additions & 0 deletions math/src/test/java/com/powsybl/math/matrix/ComplexMatrixTest.java
Original file line number Diff line number Diff line change
@@ -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.));

}
}

0 comments on commit f8ef9c6

Please sign in to comment.