Dreamlike
0%

Gams结合python的使用

發表於 更新於 分類於
文章字數： 15k 所需閱讀時間 ≈ 14 分鐘

参考资料：https://www.gams.com/latest/docs/API_PY_TUTORIAL.html

引入必要模块

1
2
3
4
5
# gams模块
from gams import *
# 常用python自带模块
import os
import sys

多Gams系统下的选择

1
2
3
4
5
# 多gams版本环境下，shell中运行py文件时可选择与default不同的gams版本运行，用以下语句进行控制
if len(sys.argv) > 1:
ws = GamsWorkspace(system_directory = sys.argv[1])
else:
ws = GamsWorkspace()

数据结构

1
2
3
4
5
6
7
8
9
10
11
plants = [ "Seattle", "San-Diego" ]
markets = [ "New-York", "Chicago", "Topeka" ]
capacity = { "Seattle": 350.0, "San-Diego": 600.0 }
demand = { "New-York": 325.0, "Chicago": 300.0, "Topeka": 275.0 }
distance = { ("Seattle", "New-York") : 2.5,
("Seattle", "Chicago") : 1.7,
("Seattle", "Topeka") : 1.8,
("San-Diego", "New-York") : 2.5,
("San-Diego", "Chicago") : 1.8,
("San-Diego", "Topeka") : 1.4
}

在python代码中，gams中的set用list表示(如plants或markets)，gams中的parameter则用dic表示(如demand或capacity)。

  • GamsDatabase 是python数据与gams进行交换的中转变量。 
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    # 创建新的GamsDatabase实例
    db = ws.add_database()

    # 往GamsDatabase中增添一维集合 'i' ，文本'canning plants'进行解释说明
    i = db.add_set("i", 1, "canning plants")
    for p in plants:
    i.add_record(p)

    # 增加定义域为 i' 的参数 'a'
    a = db.add_parameter_dc("a", [i], "capacity of plant i in cases")
    for p in plants:
    a.add_record(p).value = capacity[p]

    # 在项目空间的工作目录中，将GamsDatabase输出到名为'data.gdx' 的GDX （GAMS Data eXchange）文件
    db.export("data.gdx")
    GDX文件可在gams studio中打开并显示。同时可使用gdx-pandas将GDX转化为DataFrame结构。

运行gms文件

1
2
3
# 注意文件要有路径
t1 = ws.add_job_from_file(路径/"你的文件.gms")
t1.run()

可在运行前进行option设置。 

1
2
3
4
# 将xpress设置为默认solver
opt = ws.add_options()
opt.all_model_types = "xpress"
t1.run(opt)

输出运行日志

1
2
3
4
5
6
7
8
9
10
path=os.getcwd()# 获取当前路径
file = open(r"%s\xpress.opt"%(path), "w")# 打开一个xpress.opt设置文件
file.write("algorithm=barrier")# 在设置文件中进行算法设置
file.close()
opt.optfile = 1# 改设置令solver运行前读取opt设置
# 运行日志输出到本地文件xpress.log中
with open("xpress.log", "w") as log:
t1.run(opt, output=log)
# 如果想输出到显示屏上
t1.run(opt, output=sys.stdout)

输出运行结果

以如下小例子为例。gms代码如下 

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
Set
i /i1*i3/;

parameter
C(i)
/
i1 1
i2 2
i3 3
/;
free variable
x(i)
Z;

equations
objective
cons;

objective.. Z =e= sum(i,C(i)*x(i));
cons(i).. x(i) =g= 0;

Model MyModel /all/;

solve MyModel using lp minimizing Z;
display x.l, x.m;
pyton代码 
1
2
3
4
path=os.getcwd()
t1 = ws.add_job_from_file(r"%s\test.gms"%(path))
 for rec in t1.out_db.get_symbol("x"):
        print("x(" + rec.keys[0] +  "): level=" + str(rec.level) + " marginal=" + str(rec.marginal))
其中t1.out_db是GamsDatabase实例，也就可以输出到GDX文件中，然后用gams studio打开显示结果。

模型与数据分离(transport2.py为例)

代码中首先定义的数据函数，返回遵循gams语法规则的相关输入数据的定义。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def get_data_text():
    return '''
  Sets
       i   canning plants   / seattle, san-diego /
       j   markets          / new-york, chicago, topeka / ;

  Parameters

       a(i)  capacity of plant i in cases
         /    seattle     350
              san-diego   600  /

       b(j)  demand at market j in cases
         /    new-york    325
              chicago     300
              topeka      275  / ;

  Table d(i,j)  distance in thousands of miles
                    new-york       chicago      topeka
      seattle          2.5           1.7          1.8
      san-diego        2.5           1.8          1.4  ;

  Scalar f  freight in dollars per case per thousand miles  /90/ ; '''
接着又定义了模型函数，返回遵循gams语法规则的模型文本。其中只给出了集合及参数的声明，而未给出具体数值。 
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
def get_model_text():
    return '''
  Sets
       i   canning plants
       j   markets

  Parameters
       a(i)   capacity of plant i in cases
       b(j)   demand at market j in cases
       d(i,j) distance in thousands of miles
  Scalar f  freight in dollars per case per thousand miles;

$if not set incname $abort 'no include file name for data file provided'
$include %incname%

  Parameter c(i,j)  transport cost in thousands of dollars per case ;

            c(i,j) = f * d(i,j) / 1000 ;

  Variables
       x(i,j)  shipment quantities in cases
       z       total transportation costs in thousands of dollars ;

  Positive Variable x ;

  Equations
       cost        define objective function
       supply(i)   observe supply limit at plant i
       demand(j)   satisfy demand at market j ;

  cost ..        z  =e=  sum((i,j), c(i,j)*x(i,j)) ;

  supply(i) ..   sum(j, x(i,j))  =l=  a(i) ;

  demand(j) ..   sum(i, x(i,j))  =g=  b(j) ;

  Model transport /all/ ;

  Solve transport using lp minimizing z ;

  Display x.l, x.m ; '''
其中以下两句判断是否有名为incname的option被引入。经测试，可放置在get_model_text()的return文本中的任意位置。 
1
2
$if not set incname $abort 'no include file name for data file provided'
$include %incname%
python执行代码如下： 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
if __name__ == "__main__":
    if len(sys.argv) > 1:
        ws = GamsWorkspace(system_directory = sys.argv[1])
    else:
        ws = GamsWorkspace()

    # 将文本输出到tdata.gms文件中  注意这里的路径为ws.working_directory，最好不要自己该路径
    file = open(os.path.join(ws.working_directory, "tdata.gms"), "w")
    file.write(get_data_text())
    file.close()

    # 将字符串定义的模型添加到job中
    t2 = ws.add_job_from_string(get_model_text())
    # 添加设置
    opt = ws.add_options()
    # 定义名为“incnmae”的option(可随意选择名称，但需要与上面$if...中的一致)。“tdata”是刚输出的data_text文件，也可写为"tdata.gms"或“路径/tdata.gms”
    opt.defines["incname"] = "tdata"
    # 执行代码并输出
    t2.run(opt)
    for rec in t2.out_db["x"]:
        print("x(" + rec.key(0) + "," + rec.key(1) + "): level=" + str(rec.level) + " marginal=" + str(rec.marginal))

尝试改动，将变量定义放到get_data_text中，然后又定义了get_variable_set函数固定x的值。在get_model_text中多加一个include。在运行代码中添加了option。 

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
def get_data_text():
    return '''
  Sets
       i   canning plants   / seattle, san-diego /
       j   markets          / new-york, chicago, topeka / ;

  Parameters

       a(i)  capacity of plant i in cases
         /    seattle     350
              san-diego   600  /

       b(j)  demand at market j in cases
         /    new-york    325
              chicago     300
              topeka      275  / ;

  Table d(i,j)  distance in thousands of miles
                    new-york       chicago      topeka
      seattle          2.5           1.7          1.8
      san-diego        2.5           1.8          1.4  ;

  Scalar f  freight in dollars per case per thousand miles  /90/ ;

  Variables
       x(i,j)  shipment quantities in cases
       z       total transportation costs in thousands of dollars ;

  Positive Variable x ;

   '''

def get_variable_set():
    return '''
    x.fx('seattle','new-york')=25;
    '''


def get_model_text():
    return '''
$if not set incname1 $abort 'no include file name for data file provided'
$if not set incname2 $abort 'no include file name for data file provided'
$include %incname1%
$include %incname2%
  Sets
       i   canning plants
       j   markets

  Parameters
       a(i)   capacity of plant i in cases
       b(j)   demand at market j in cases
       d(i,j) distance in thousands of miles
  Scalar f  freight in dollars per case per thousand miles;

  Parameter c(i,j)  transport cost in thousands of dollars per case ;

            c(i,j) = f * d(i,j) / 1000 ;



  Equations
       cost        define objective function
       supply(i)   observe supply limit at plant i
       demand(j)   satisfy demand at market j ;

  cost ..        z  =e=  sum((i,j), c(i,j)*x(i,j)) ;

  supply(i) ..   sum(j, x(i,j))  =l=  a(i) ;

  demand(j) ..   sum(i, x(i,j))  =g=  b(j) ;

  Model transport /all/ ;

  Solve transport using lp minimizing z ;

  Display x.l, x.m ; '''

 if __name__ == "__main__":
    if len(sys.argv) > 1:
        ws = GamsWorkspace(system_directory = sys.argv[1])
    else:
        ws = GamsWorkspace()

    file = open(os.path.join(ws.working_directory, "tdata1.gms"), "w")
    file.write(get_data_text())
    file.close()

    file = open(os.path.join(ws.working_directory, "tdata2.gms"), "w")
    file.write(get_variable_set())
    file.close()

    t2 = ws.add_job_from_string(get_model_text())
    opt = ws.add_options()
    opt.defines["incname1"] = "tdata1.gms"
    opt.defines["incname2"] = "tdata2.gms"
    t2.run(opt)
    for rec in t2.out_db["x"]:
        print("x(" + rec.key(0) + "," + rec.key(1) + "): level=" + str(rec.level) + " marginal=" + str(rec.marginal))

利用GamsModifier修改模型实例的parameter和variable

首先需要了解 GamsCheckpoint ，其作用是捕捉GamsJob的状态。详细解释如下： 

1
2
3
4
5
6
# 定义一个checkpoint
cp = ws.add_checkpoint()
# 利用string定义一个GamsJob t5 但string中不含有solve语句
t5 = ws.add_job_from_string(get_model_text())
# 利用checkpoint捕捉及控制GamsJob
t5.run(checkpoint=cp)
其中get_model_text()具体如下，注意没有像前面的代码中一样包含“ Solve transport using lp minimizing z ;”的语句。 
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
def get_model_text():
    return '''
  Sets
       i   canning plants   / seattle, san-diego /
       j   markets          / new-york, chicago, topeka / ;

  Parameters

       a(i)  capacity of plant i in cases
         /    seattle     350
              san-diego   600  /

       b(j)  demand at market j in cases
         /    new-york    325
              chicago     300
              topeka      275  / ;

  Table d(i,j)  distance in thousands of miles
                    new-york       chicago      topeka
      seattle          2.5           1.7          1.8
      san-diego        2.5           1.8          1.4  ;

  Scalar f      freight in dollars per case per thousand miles  /90/ ;
  Scalar bmult  demand multiplier /1/;

  Parameter c(i,j)  transport cost in thousands of dollars per case ;

            c(i,j) = f * d(i,j) / 1000 ;

  Variables
       x(i,j)  shipment quantities in cases
       z       total transportation costs in thousands of dollars ;

  Positive Variable x ;

  Equations
       cost        define objective function
       supply(i)   observe supply limit at plant i
       demand(j)   satisfy demand at market j ;

  cost ..        z  =e=  sum((i,j), c(i,j)*x(i,j)) ;

  supply(i) ..   sum(j, x(i,j))  =l=  a(i) ;

  demand(j) ..   sum(i, x(i,j))  =g=  bmult*b(j) ;

  Model transport /all/ ;
  Scalar ms 'model status', ss 'solve status'; '''
利用checkpoint可以实现变换输入，进行多次Job运算，得到相应的结果。 
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
if __name__ == "__main__":
    if len(sys.argv) > 1:
        ws = GamsWorkspace(system_directory = sys.argv[1])
    else:
        ws = GamsWorkspace()

    cp = ws.add_checkpoint()

    # initialize a GAMSCheckpoint by running a GAMSJob
    t5 = ws.add_job_from_string(get_model_text())
    t5.run(checkpoint=cp)

    # 定义了bmultlist的不同取值
    bmultlist = [ 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3 ]

    # create a new GAMSJob that is initialized from the GAMSCheckpoint
    for b in bmultlist:
        # 针对不同的bmult取值生成具有不同输入的Job t5，并添加solve 语句。checkpoint cp 也作为参数输入，个人理解是新输入的语句+cp拼合成完整包含不同参数和求解语句的模型。
        t5 = ws.add_job_from_string("bmult=" + str(b) + "; solve transport min z use lp; ms=transport.modelstat; ss=transport.solvestat;", cp)
        # 执行Job并输出结果，结果储存在GamsDatabase类型变量t5.out_db中。
        t5.run()
        print("Scenario bmult=" + str(b) + ":")
        print("  Modelstatus: " + str(t5.out_db["ms"].find_record().value))
        print("  Solvestatus: " + str(t5.out_db["ss"].find_record().value))
        print("  Obj: " + str(t5.out_db["z"].find_record().level))
接着介绍如何使用checkpoint创建GamsModelInstance。个人感觉GamsModelInstance的作用跟Job是一样的，都是可以用来进行模型运算，输出结果。同时可以实现变动参数和变动变量下的模型输出。 
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
if __name__ == "__main__":
    if len(sys.argv) > 1:
        ws = GamsWorkspace(system_directory = sys.argv[1])
    else:
        ws = GamsWorkspace()
    # 生成checkpoint
    cp = ws.add_checkpoint()

    # initialize a GAMSCheckpoint by running a GAMSJob
    t7 = ws.add_job_from_string(get_model_text())# 生成包含模型部分的job，但不包含solve语句
    t7.run(checkpoint=cp)# 运行job，将状态传递给checkpoint，大约意味着checkpoint获得了模型信息，包括数据和公式

    # 下面部分是实现每次运行都可以有不同的参数设置
    # create a GAMSModelInstance and solve it multiple times with different scalar bmult
    mi = cp.add_modelinstance()# 生成模型实例
    bmult = mi.sync_db.add_parameter("bmult", 0, "demand multiplier")# 向改实例中添加参数，函数参数2:0表示维度dim为0 可输出到gdx中查看，此时bmult还是空值，结果为下图1
    opt = ws.add_options()# 向工作空间中添加option，这个是公用的
    opt.all_model_types = "cplexd"# 设定求解器
    # opt.all_model_types = "gurobi" 文档中给出的另一个实例

    # instantiate the GAMSModelInstance and pass a model definition and GAMSModifier to declare bmult mutable
    mi.instantiate("transport use lp min z", GamsModifier(bmult), opt)# 参数1增添了solve语句；参数2表示参数bmult在运行时是可变的，而其他参数等式不变；参数3表示使用上面给的设置；此时输出到gdx可见模型参数已经随opt传入，bmult依旧是空值，具体见图2

    bmult.add_record().value = 1.0# 给可变参数bmult添加初始值 此时输出gdx可见bmult已经有数值，见图3
    bmultlist = [ 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3 ]

    print("Scenario bmult=" + str(mi.sync_db['bmult'].find_record().value))# 可显示当前场景下bmult的数值 下面语句一样的效果
    print("Scenario bmult=" + str(mi.sync_db.get_parameter('bmult').find_record().value))# 可显示当前场景下bmult的数值


    for b in bmultlist:
        bmult.first_record().value = b# 求解前给bmult赋不同值
        mi.solve()# 模型求解
        print("Scenario bmult=" + str(b) + ":")
        print("  Modelstatus: " + str(mi.model_status))
        print("  Solvestatus: " + str(mi.solver_status))
        print("  Obj: " + str(mi.sync_db.get_variable("z").find_record().level))

    #下面部分是实现运行时有不同的变量设置
    # create a GAMSModelInstance and solve it with single links in the network blocked
    mi = cp.add_modelinstance()# 生成实例
    x = mi.sync_db.add_variable("x", 2, VarType.Positive)# 向实例中添加变量，函数参数2：“2”表示维度是2，函数参数3：VarType.Positive表示变量为正实数
    xup = mi.sync_db.add_parameter("xup", 2, "upper bound on x")
    # 此时输出到gdx，可见变量x和参数xup都还没哟具体数值。见图4

    # instantiate the GAMSModelInstance and pass a model definition and GAMSModifier to declare upper bound of X mutable
    mi.instantiate("transport use lp min z", GamsModifier(x, UpdateAction.Upper, xup))# 生成实例。GamsModifier用来制定x是的某些属性是可更改的，参数UpdateAction.Upper即表示改属性（具体表示x的上界upper bound可变）。文档给出三种可变属性：
lower bound, upper bound or level。UpdateAction类的详细信息科参见文档，包含Upper、Lower 、Fixed、Primal和Dual等五种选择。参数xup为对应的属性更换的数值。虽然没有像上面参数的例子添加opt参数，但输出的gdx中已经包含模型的集合等数据。见图5。
    mi.solve()

    for i in t7.out_db["i"]:
        for j in t7.out_db["j"]:
            xup.clear()# 清除xup设置
            xup.add_record((i.key(0),j.key(0))).value = 0 # 为某一link设置上界为零，见图6
            mi.solve()
            print("Scenario link blocked: " + i.key(0)  + " - " + j.key(0))
            print("  Modelstatus: " + str(mi.model_status))
            print("  Solvestatus: " + str(mi.solver_status))
            print("  Obj: " + str(mi.sync_db["z"].find_record().level))

    # 尝试使用。目标是固定某link必须使用，经测试成功。
    mi = cp.add_modelinstance()
    x = mi.sync_db.add_variable("x", 2, VarType.Positive)
    xlevel = mi.sync_db.add_parameter("xlevel", 2, "fixed value of x")
    mi.instantiate("transport use lp min z", GamsModifier(x, UpdateAction.Fixed, xlevel))# 令某link的x固定为xlevel=1
    for i in t7.out_db["i"]:_config.yml
        for j in t7.out_db["j"]:
            xlevel.clear()
            xlevel.add_record((i.key(0), j.key(0))).value = 1
            # mi.solve()
            path = os.getcwd()
            mi.sync_db.export(r"%s\d.gdx" % (path))
            mi.solve()
            print("Scenario link forced: " + i.key(0) + " - " + j.key(0))
            print("  Modelstatus: " + str(mi.model_status))
            print("  Solvestatus: " + str(mi.solver_status))
            print("  Obj: " + str(mi.sync_db["z"].find_record().level))
图1： 图2： 图3: 图4： 图5： 图6：

  • 有时需要从gms文件加载Job，但如果不特别设置，ws的working_directory是一个临时文件夹，可能找不到gms文件。在 https://forum.gamsworld.org/viewtopic.php?t=10228 给出了详细的解释。解决方案是： 

    1
    ws = GamsWorkspace(working_directory=".")

  • 有时需要修改多个参数，这时就需要修改instantiate函数的第二个参数，改单一GamsModifier为跑list 

    1
    mi.instantiate("transport use lp min z", [GamsModifier(g1),GamsModifier(g2)...,GamsModifier(gn)], opt)
    参考说明文档:https://www.gams.com/latest/docs/apis/java/classcom_1_1gams_1_1api_1_1GAMSModelInstance.html#af6afd2f81e9e3ca7ee4656ee467afe49

  • 有时需要修改一个parameters，也就是一个list的初始值。这时比较麻烦，因为该参数需要set作为add函数的第二个输入变量，然而这个set无法一开始就在sync_db中定义。于是需要一个迂回的操作，具体如下：

1
2
3
4
5
6
7
8
9
10
11
12
# 得到集合 用作add_parameter的第二个参数 此时record locked无法增加修改
Si=mi.sync_db.add_set("Si",1,"集合")
S0 = mi.sync_db.add_parameter_dc("S0", [Si], "需要修改的parameters")
# set 不允许作为GamsModifier。instantiate后就可以修改record
mi.instantiate("problem using mip max Z", [其他的GamsModifier...,GamsModifier(S0)])

# 为了能修改S0 即使gms文件中已经有Si的定义 但在这里还需要定义一遍
for s in 给出的Set:#给出的Set  就是一个list 如['i1','i2',...,'in']
    Si.add_record(s)
# 给S0赋初始值
for s in S0的索引:
    S0.add_record(s).value = 具体要赋予的初始值数值