gamspy-examples/models/syscomp/syscomp.py at master · GAMS-dev/gamspy-examples · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
"""
## 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 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()

    # 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()