Solving complex linear algebraic systems of equations

syscomp.py

"""
## LICENSETYPE: Demo
## MODELTYPE: LP


Solving complex linear algebraic systems of equations

References
Kalvelagen, E., (2002) Solving systems of linear equations with GAMS.
http://www.gams.com/~erwin/lineq.pdf
"""

from __future__ import annotations

import os

import numpy as np
import pandas as pd
from gamspy import (
    Alias,
    Container,
    Equation,
    Model,
    Parameter,
    Set,
    Sum,
    Variable,
)


def data_records():
    cols = ["i1", "i2", "rhs"]
    idxs = [("real", "i1"), ("real", "i2"), ("imag", "i1"), ("imag", "i2")]

    data = np.array([[30, 20, 14], [15, 8, 11], [10, -15, 5], [0, -4, -7]])

    idxs = pd.MultiIndex.from_tuples(idxs, names=["Index1", "Index2"])
    data = pd.DataFrame(data, columns=cols, index=idxs)
    data.reset_index(inplace=True)
    melted_data = data.melt(
        id_vars=["Index1", "Index2"], value_vars=["i1", "i2", "rhs"]
    )
    return melted_data


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
    )

    # SET #
    i = Set(m, name="i", records=["i1", "i2"])

    # ALIAS #
    j = Alias(m, name="j", alias_with=i)

    # PARAMETER #
    data = Parameter(
        m, name="data", domain=["*", "*", "*"], records=data_records()
    )

    # VARIABLES #
    rx = Variable(
        m, name="rx", domain=i, description="real part of the solution"
    )
    ix = Variable(
        m, name="ix", domain=i, description="imaginary part of the solution"
    )

    # EQUATIONS #
    real = Equation(
        m,
        name="real",
        type="regular",
        domain=i,
        description="real part of the system",
    )
    imag = Equation(
        m,
        name="imag",
        type="regular",
        domain=i,
        description="imaginary part of the system",
    )

    real[i] = (
        Sum(j, data["real", i, j] * rx[j] - data["imag", i, j] * ix[j])
        == data["real", i, "rhs"]
    )
    imag[i] = (
        Sum(j, data["imag", i, j] * rx[j] + data["real", i, j] * ix[j])
        == data["imag", i, "rhs"]
    )

    syscomp = Model(
        m,
        name="syscomp",
        equations=m.getEquations(),
        problem="lp",
        sense="FEASIBILITY",
    )
    syscomp.solve()

    # REPORTING PARAMETER
    rep = Parameter(m, name="rep", domain=["*", i])

    rep["rx", i] = rx.l[i]
    rep["ix", i] = ix.l[i]

    print(
        "Objective Function Value:  ", round(syscomp.objective_value, 4), "\n"
    )
    print("Solution Summary:\n", rep.pivot().round(3))

    # End of SysComp


if __name__ == "__main__":
    main()