# A General MINLP Model for the Multiway Valve Channel-Limited Location-Allocation Problem

Authors: Boyu Zhu; Jun Zhou; Guangchuan Liang; Xuan Zhou; Liulin Zhou
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

As an important aspect of oil and gas field construction, oil and gas field surface engineering often affects the efficiency and safety of oil and gas production, and it requires a large investment. In the past, much work has been done on layout optimization for gathering pipeline networks. However, few studies have considered multiway valves in the optimization of pipeline networks. Compared to traditional metering processes, a process using multiway valves can reduce construction and operation costs and enable automation of the well-selection operation. In this paper, an MINLP model is established in which the number of multiway valves and their numbers of channels are considered special constraints and the number of multiway valves and the associations between wells and multiway valves are treated as optimization variables. A specific heuristic algorithm for solving this problem is also proposed in this paper. We consider the coordinates of real-world wells and of randomly generated well locations as different examples to analyze the performance of this algorithm. These examples demonstrate that the algorithm can initially adjust the network through single-step iteration, while double-step iteration is efficient when all channels of a multiway valve have been associated with pipelines, and multistep iteration can help the objective function escape from local optima. Finally, numerical analysis results prove that the proposed algorithm can be used to efficiently solve the problem of interest and exhibits stable convergence.


---

## Body

## 1. Introduction

In recent years, as the production scale of oil and gas fields has expanded, the construction and labor costs incurred for the measurement of oil and gas production per well have also increased. In addition, the practical implementation of conventional metering methods is subject to certain limitations, especially in areas such as deserts and deep-sea regions. There are three main metering methods for oil and gas fields: single-well metering, metering by manual rotation, and automatic rotation metering. The single-well metering process requires the installation of metering equipment for each well, as shown in Figure1(a). The initial investment is relatively expensive, especially when the well production rate is below the initial estimate. In the manual rotation measurement process, as shown in Figure 1(b), each valve that is to be opened or closed needs to be operated manually. Although it is not necessary to install metering devices at each well site, metering a well will require the manual operation of more than a dozen valves, and after metering, these valves will need to be operated again to discharge the liquid in the metering separator, thus requiring more operation time and incurring a higher cost. In addition, the metering device occupies a large area, and such devices are not suitable for offshore platforms. Regarding the third method, retrofitting an actuator onto the original manual valve of an existing manifold will greatly increase the construction cost. Therefore, it is necessary to find a device that can enable automatic well-selection metering while incurring lower construction and operation costs.Figure 1
Comparison of metering processes: (a) single-well metering; (b) metering by manual rotation; (c) automatic rotation metering.
(a)
(b)
(c)At present, many oil and gas field production enterprises use multiway valves (MVs) to achieve automatic rotation metering. The MV-based metering process can reduce the amount of manual operation required in remote areas, which is greatly beneficial for reducing operation costs. The metering process with an MV is shown in Figure1(c). The blue dashed arrow indicates the flow that needs to be metered at present, and the black arrows indicate the flows from production wells. When switching to metering another well, the control system requires only remote operation, thus leading to considerable savings in operation cost. Compared to the traditional methods, the MV-based metering process reduces the numbers of valves, gauges, and actuators that need to be installed on the pipelines. Each device also occupies a smaller space, making these devices more suitable for a wider variety of oil and gas fields, especially for coalbed methane mining, shale gas, and other unconventional oil and gas resources. From here on, when the term “valve” is used in this paper, it refers to an MV.The number of channels of an MV is a characteristic that is defined by the manufacturer and determines how many pipelines can access the MV. If an MV has 16 channels, this means that at most 16 wells can be connected to the MV. Manufacturers make MVs with various channel counts, such as 12, 14, and 16. Depending on the conditions in the field or other requirements such as consistency principle, the designer needs to select MVs with an appropriate number of channels to meet the requirements of the project. Therefore, during the design process, the number of wells connected to an MV is constrained by the number of channels of the MV, which depends on the designer’s selection. Depending on the constraint relaxation degree, there are three possible ways to specify MV-related constraints:(1)
The number of MVs is unrestricted(2)
Allowable ranges are given for the number of MVs and their channel counts(3)
The channel counts of the MVs and the number of MVs are fixed in advanceFor a given target oil and gas reservoir, the traditional solution is to invite experts to design feasible solutions based on real scenarios, but this approach does not necessarily guarantee the most economical cost. By contrast, optimization methods can help a designer to obtain the best solution. However, in the past, researchers have rarely considered constraints related to the channel counts of the MVs and the number of MVs, although such constraints are important to ensure that the optimal solution satisfies the required service conditions. In this paper, the problem of interest is formulated as a special P-median problem, and the constraints mentioned above are integrated into the mathematical model. By treating the associations between wells and MVs and the number of MVs as variables and the minimum cost as the objective function, a general mixed integer nonlinear programming (MINLP) model is established. Based on the connection and foresight matrices of the production wells and MVs, a foresight matrix search algorithm (FMSA) is proposed in this paper to solve the MINLP problem. In this algorithm, the nonlinear part of the problem is treated as a subproblem, and the optimization results are obtained through iteration. By examples, the convergence of the algorithm can be ensured, as can the quality and efficiency of the solution. In addition, it is proven that, with this algorithm, the solution process can be completed in a reasonable time even if the problem is large in scale.This paper is divided into six sections: In Section1, we have introduced the problem motivating this paper and what we have done to address it. In Section 2, we review the related research contributions that have been reported in the past. In Section 3, the detailed mathematical model is presented, and the special constraints are formulated. Subsequently, the details of solving the mathematical model using the proposed algorithm are introduced in Section 4. The results are analyzed in Section 5. Finally, conclusions and discussions are presented in Section 6. The data used in this paper can be found in Tables 1–5.Table 1
Coordinates of 192 wells from a gas field.
x (m)
y (m)
Number
x (m)
y (m)
Number
x (m)
y (m)
Number
7430
2770
W1
4176
5039
W65
5570
3473
W129
7177
2698
W2
3975
4957
W66
5801
3634
W130
7193
3084
W3
4068
5358
W67
4985
3832
W131
7112
3341
W4
4429
5300
W68
7389
5655
W132
7572
2515
W5
4722
5370
W69
8522
5800
W133
8827
4759
W6
4904
5150
W70
8103
5715
W134
9009
5366
W7
4273
5543
W71
7542
5441
W135
9587
5795
W8
4558
5660
W72
7186
5593
W136
9458
5509
W9
5138
5270
W73
6566
5406
W137
9135
5732
W10
4969
5575
W74
6475
4116
W138
8935
5562
W11
5231
5060
W75
6024
4415
W139
8746
5209
W12
6030
5549
W76
6221
4302
W140
8774
5516
W13
5468
4806
W77
6465
4443
W141
8565
5454
W14
6472
5940
W78
6295
4625
W142
8261
4764
W15
7846
5656
W79
6120
4893
W143
8253
4440
W16
6208
5351
W80
6374
5033
W144
8396
4255
W17
5788
5334
W81
6650
5172
W145
8479
4649
W18
5557
5251
W82
6523
4829
W146
8713
4407
W19
5749
4979
W83
6711
4567
W147
9132
4226
W20
5994
5116
W84
5997
3457
W148
8857
4207
W21
6756
4798
W85
5990
3836
W149
8479
4998
W22
5723
4398
W86
6292
3902
W150
722
1355
W23
5573
4520
W87
6206
3526
W151
1779
1152
W24
5965
2209
W88
6492
3690
W152
1497
1374
W25
6276
2157
W89
6760
3709
W153
1971
1398
W26
6395
2580
W90
6593
3934
W154
1952
952
W27
6095
2413
W91
6250
2881
W155
1617
989
W28
5827
2630
W92
5532
2840
W156
1671
763
W29
5791
2352
W93
5819
2923
W157
1352
1125
W30
5384
2478
W94
5774
3200
W158
890
1183
W31
5587
2195
W95
6151
3119
W159
535
738
W32
5421
1986
W96
6405
3284
W160
264
29
W33
5718
2003
W97
6638
3422
W161
140
460
W34
6953
1171
W98
6812
3137
W162
459
285
W35
6758
1383
W99
6962
2906
W163
284
2160
W36
6489
1199
W100
6693
2722
W164
1805
2128
W37
6641
978
W101
6525
2993
W165
1466
1868
W38
6971
904
W102
7089
4415
W166
4185
4644
W39
7321
1026
W103
7254
4146
W167
4266
4411
W40
7397
1444
W104
7067
4733
W168
3996
4274
W41
7178
1236
W105
7036
5093
W169
4167
4079
W42
7074
1509
W106
7340
5187
W170
4047
3589
W43
6615
1674
W107
7832
5149
W171
4213
3355
W44
6959
2201
W108
8062
5110
W172
4280
3785
W45
6836
2432
W109
7609
4981
W173
4439
4188
W46
6529
2291
W110
7273
4862
W174
4642
3905
W47
6646
2086
W111
7835
4858
W175
4506
3500
W48
6450
1842
W112
7626
4730
W176
4745
3665
W49
6228
1732
W113
7819
4479
W177
4902
3436
W50
7117
1848
W114
7883
4211
W178
4624
3289
W51
7220
2425
W115
8024
4616
W179
4849
4060
W52
6837
1771
W116
7521
4356
W180
4417
4898
W53
6076
1312
W117
7333
4573
W181
4553
4667
W54
6304
1463
W118
7031
3908
W182
4614
4356
W55
5601
1025
W119
7409
3886
W183
4767
4755
W56
5263
808
W120
7728
4028
W184
4910
4514
W57
5429
1272
W121
7150
3733
W185
4671
5018
W58
5809
1246
W122
7370
3515
W186
5016
4291
W59
5721
1462
W123
7599
3580
W187
5021
4926
W60
5238
1591
W124
7878
3778
W188
5243
4666
W61
5533
1698
W125
7735
3369
W189
5220
4048
W62
5887
1725
W126
7497
3253
W190
5352
4367
W63
5219
3276
W127
7883
3225
W191
5465
4197
W64
5441
3674
W128
7696
2928
W192Table 2
Coordinates of 35 wells from a gas field.
x (m)
y (m)
Number
x (m)
y (m)
Number
x (m)
y (m)
Number
264
29
W1
4213
3355
W13
4767
4755
W25
140
460
W2
4280
3785
W14
4910
4514
W26
459
285
W3
4439
4188
W15
4671
5018
W27
284
2160
W4
4642
3905
W16
5016
4291
W28
1805
2128
W5
4506
3500
W17
5021
4926
W29
1466
1868
W6
4745
3665
W18
5243
4666
W30
1793
2404
W7
4902
3436
W19
5220
4048
W31
4185
4644
W8
4624
3289
W20
5352
4367
W32
4266
4411
W9
4849
4060
W21
5465
4197
W33
3996
4274
W10
4417
4898
W22
4176
5039
W34
4167
4079
W11
4553
4667
W23
3975
4957
W35
4047
3589
W12
4614
4356
W24Table 3
Coordinates of 40 wells randomly generated by MATLAB.
x (m)
y (m)
Number
x (m)
y (m)
Number
x (m)
y (m)
Number
855.158
2624.822
W1
1365.531
7212.275
W15
9786.806
7126.945
W29
8010.146
292.2028
W2
1067.619
6537.573
W16
5004.716
4710.884
W30
9288.541
7303.309
W3
4941.739
7790.517
W17
596.1887
6819.719
W31
4886.09
5785.251
W4
7150.371
9037.206
W18
424.3114
714.4546
W32
2372.836
4588.488
W5
8909.225
3341.631
W19
5216.498
967.3003
W33
9630.885
5468.057
W6
6987.458
1978.098
W20
8181.486
8175.471
W34
5211.358
2315.944
W7
305.4095
7440.743
W21
7224.396
1498.654
W35
4888.977
6240.601
W8
5000.224
4799.221
W22
6596.053
5185.949
W36
6791.355
3955.152
W9
9047.222
6098.666
W23
9729.746
6489.915
W37
3674.366
9879.82
W10
6176.664
8594.423
W24
8003.306
4537.977
W38
377.3887
8851.68
W11
8054.894
5767.215
W25
4323.915
8253.138
W39
9132.868
7961.839
W12
1829.225
2399.32
W26
834.6981
1331.71
W40
987.1228
2618.712
W13
8865.119
286.7415
W27
3353.568
6797.28
W14
4899.014
1679.271
W28Table 4
Coordinates of 80 wells randomly generated by MATLAB.
x (m)
y (m)
Number
x (m)
y (m)
Number
x (m)
y (m)
Number
1269.868
9133.759
W1
6554.779
1711.867
W28
5383.424
9961.347
W55
1576.131
9705.928
W2
7060.461
318.328
W29
8173.032
8686.947
W56
4217.613
9157.355
W3
6550.980
1626.117
W30
8692.922
5797.046
W57
971.318
8234.578
W4
5852.678
2238.119
W31
8530.311
6220.551
W58
3516.595
8308.286
W5
5472.155
1386.244
W32
9027.161
9447.872
W59
1621.823
7942.845
W6
4505.416
838.214
W33
2784.982
5468.815
W60
2289.770
9133.374
W7
5498.602
1449.548
W34
1868.726
4897.644
W61
1523.780
8258.170
W8
4387.444
3815.585
W35
2760.251
6797.027
W62
1066.528
9618.981
W9
4455.862
6463.130
W36
1189.977
4983.641
W63
46.342
7749.105
W10
5059.571
6990.767
W37
119.021
3371.226
W64
2598.704
8000.685
W11
6160.447
4732.888
W38
3112.150
5285.331
W65
4314.138
9106.476
W12
5852.641
5497.236
W39
1656.487
6019.819
W66
3377.194
9000.538
W13
3804.458
5678.216
W40
2629.713
6540.791
W67
9571.669
4853.756
W14
5688.237
4693.906
W41
781.755
4426.783
W68
8002.805
1418.863
W15
3509.524
5132.495
W42
844.358
3997.826
W69
7431.325
3922.270
W16
4908.641
4892.526
W43
2769.230
461.714
W70
6948.286
3170.995
W17
8147.237
9057.919
W44
1492.940
2575.083
W71
9502.220
344.461
W18
9575.068
9648.885
W45
1965.953
2510.839
W72
9597.440
3403.857
W19
7922.073
9594.924
W46
758.543
539.501
W73
7512.671
2550.951
W20
8491.293
9339.932
W47
1818.470
2638.029
W74
8407.173
2542.822
W21
6787.352
7577.401
W48
1455.390
1360.686
W75
8142.848
2435.250
W22
7655.168
7951.999
W49
4018.080
759.667
W76
9292.636
3499.838
W23
7093.648
7546.867
W50
2399.162
1233.189
W77
9171.937
2858.390
W24
8909.033
9592.914
W51
1839.078
2399.525
W78
9340.107
1299.062
W25
7572.002
7537.291
W52
4172.671
496.544
W79
6323.592
975.404
W26
5307.976
7791.672
W53
3692.468
1112.028
W80
6557.407
357.117
W27
6892.145
7481.516
W54Table 5
Coordinates of 191 wells randomly generated by MATLAB.
x (m)
y (m)
Number
x (m)
y (m)
Number
x (m)
y (m)
Number
5438.859
7210.466
W1
6951.63
4991.16
W65
2547.902
2240.4
W129
3658.162
7635.046
W2
377.3887
8851.68
W66
1886.62
2874.982
W130
3934.564
6714.311
W3
1365.531
7212.275
W67
2702.704
1970.538
W131
5522.913
6298.834
W4
1067.619
6537.573
W68
2702.943
2084.614
W132
4670.682
6481.984
W5
305.4095
7440.743
W69
2955.338
3329.363
W133
4241.668
5078.583
W6
596.1887
6819.719
W70
8443.088
1947.643
W134
4886.09
5785.251
W7
154.8713
9840.637
W71
8010.146
292.2028
W135
5000.224
4799.221
W8
526.77
7378.581
W72
6987.458
1978.098
W136
5004.716
4710.884
W9
326.0082
5611.998
W73
8865.119
286.7415
W137
5268.758
4167.995
W10
1564.05
8555.228
W74
7224.396
1498.654
W138
4896.876
3394.934
W11
911.1346
5762.094
W75
9990.804
1711.211
W139
5985.237
4709.243
W12
444.5409
7549.333
W76
8177.606
2607.28
W140
4400.851
5271.427
W13
5.223754
8654.386
W77
6951.405
679.9277
W141
4624.492
4243.49
W14
834.8281
6259.598
W78
9159.912
11.51057
W142
5822.492
5407.393
W15
319.9102
6147.134
W79
8699.41
2647.79
W143
4045.8
4483.729
W16
426.5241
6351.979
W80
6962.663
938.2003
W144
1671.684
1062.163
W17
252.2818
8422.066
W81
7688.543
1672.535
W145
1781.325
1280.144
W18
3507.271
9390.016
W82
8013.476
2278.429
W146
1614.847
1787.662
W19
3111.023
9233.796
W83
8448.557
2094.051
W147
2362.306
1193.962
W20
3674.366
9879.82
W84
9479.331
820.7121
W148
2186.766
1057.983
W21
4323.915
8253.138
W85
7378.417
634.045
W149
1096.975
635.9137
W22
5478.709
9427.37
W86
4709.233
2304.882
W150
1920.283
1388.742
W23
4177.441
9830.525
W87
4302.074
1848.163
W151
1117.057
1362.925
W24
4607.259
9816.38
W88
4388.7
1111.192
W152
1476.082
549.7415
W25
4574.244
8753.716
W89
5211.358
2315.944
W153
169.8294
1208.596
W26
5180.521
9436.226
W90
4899.014
1679.271
W154
1057.094
1420.411
W27
5224.953
9937.046
W91
3724.097
1981.184
W155
1339.313
308.8955
W28
4386.45
8335.006
W92
4252.593
3127.189
W156
2276.643
4356.987
W29
5144.235
8842.81
W93
4228.857
942.2934
W157
2580.647
4087.198
W30
5827.91
8153.972
W94
4734.86
1527.212
W158
2372.836
4588.488
W31
6125.665
9899.502
W95
3180.741
1192.145
W159
1733.886
3909.378
W32
4980.943
9008.525
W96
3477.127
1499.973
W160
2919.841
4316.512
W33
5746.612
8451.782
W97
3773.955
2160.189
W161
2691.194
4228.356
W34
5216.498
967.3003
W98
3624.115
495.3258
W162
1904.333
3689.165
W35
4820.221
1206.116
W99
4895.7
1925.104
W163
1909.237
4282.53
W36
5895.075
2261.877
W100
4170.29
2059.755
W164
2427.854
4424.023
W37
5840.693
1077.69
W101
4888.977
6240.601
W165
3308.579
4243.095
W38
5943.563
225.1259
W102
4941.739
7790.517
W166
1998.628
4069.548
W39
6385.308
336.0384
W103
3846.191
5829.864
W167
1897.104
4950.058
W40
6753.321
67.15314
W104
5308.643
6544.457
W168
2818.669
5385.967
W41
4713.572
357.6273
W105
4076.192
8199.812
W169
1239.323
4903.573
W42
4423.054
195.7762
W106
4609.164
7701.597
W170
1776.025
3985.895
W43
5880.261
1547.523
W107
3411.246
6073.892
W171
855.158
2624.822
W44
5340.641
899.5068
W108
4257.288
6444.428
W172
987.1228
2618.712
W45
5737.098
520.7789
W109
4587.255
6619.448
W173
604.7118
3992.578
W46
6224.751
5870.447
W110
4794.632
6393.17
W174
688.061
3195.997
W47
6028.431
7112.158
W111
5447.161
6473.115
W175
908.2329
2664.715
W48
6596.053
5185.949
W112
9132.868
7961.839
W176
1536.567
2810.053
W49
6568.599
6279.734
W113
7150.371
9037.206
W177
1230.837
2054.942
W50
6663.389
5391.265
W114
6176.664
8594.423
W178
1465.149
1890.722
W51
6981.055
6665.279
W115
8181.486
8175.471
W179
6443.181
3786.094
W52
6959.493
6998.878
W116
8313.797
8033.644
W180
5948.961
2622.117
W53
6833.632
5465.931
W117
7183.589
9686.493
W181
6791.355
3955.152
W54
6476.176
6790.168
W118
6377.091
9576.939
W182
6447.645
3762.722
W55
6278.964
7719.804
W119
6678.327
8443.922
W183
6170.909
2652.809
W56
6834.159
7040.474
W120
6357.867
9451.741
W184
7302.488
3438.77
W57
7386.403
5859.87
W121
8419.292
8329.168
W185
5313.339
3251.457
W58
6609.446
7297.519
W122
8085.141
7550.771
W186
6021.705
3867.712
W59
7690.291
5814.465
W123
7487.057
8255.838
W187
6073.039
4501.377
W60
5649.796
6403.118
W124
2259.218
1707.08
W188
6620.096
4161.586
W61
2077.423
3012.463
W125
2217.467
1174.177
W189
5860.921
2621.453
W62
2966.759
3187.783
W126
424.3114
714.4546
W190
6877.961
3592.282
W63
1829.225
2399.32
W127
834.6981
1331.71
W191
6786.523
4951.77
W64
2518.061
2904.407
W128


## 2. Literature Review

The complexity and large scale of the corresponding mathematical models make MINLP problems a popular and difficult topic to address for scholars around the world. The methods for solving MINLP problems can be generally divided into deterministic methods and heuristic methods. Methods of the former type can yield globally optimal solutions for small-scale problems but are often inefficient for large-scale problems, whereas the quality of the solutions obtained with the heuristic methods that have been proposed to date usually cannot be guaranteed, and these methods cannot be used to solve large-scale problems within a limited time.The problem of interest here is different from the well group partitioning (WGP) problem, which usually considers factors such as the gathering radius and pipeline diameter. The purpose of the conventional WGP model is to determine the optimal mapping between wells and stations. In other words, subject to well constraints and other constraints, the wells are divided among different gathering stations to save construction costs. This problem is a kind of location-allocation problem. In a discrete location-allocation problem, P facilities are selected from a set of available candidate sites such that a given objective function is minimized. This problem is specifically referred to as the P-median problem. However, there are some differences between the model established in this paper and the general P-median model. In the problem considered here, the alternative facility (MV) locations are the same as the well locations, meaning thatm=n, whereas in the conventional P-median problem, m is often less than n. Moreover, the optimization target of the conventional P-median problem is usually to minimize the total distance from facilities to customers, while in this paper, the target is to minimize the cost of all equipment, including pipelines and devices, and the model is constrained by the special constraints of MVs. In addition, special constraints on the numbers of channels and facilities (MVs) are imposed in the model considered here. All of these factors increase the scale of the optimization problem, making the problem more complicated and difficult to solve. In [1, 2], the authors have proven that the P-median problem is an NP-hard problem when the number of facilities is uncertain. That is, the MINLP problem solved in this paper is a special case of the location-allocation problem. In this problem, the number of MVs m and the associations between wells and valves are treated as optimization variables. The associations between pipelines and various facilities are precisely the elements of interest in oil and gas field surface engineering, and reasonable optimization of the pipeline network structures between facilities can effectively reduce the engineering and construction costs. There are many possible pipeline network structures, including star, tree, and ring structures. With the development of methods for exploiting unconventional oil and gas resources, tree-tree network structures [3, 4] have been used in marginal gas reservoirs for resources such as coalbed gas and shale gas. Such a structure can considerably reduce costs by virtue of the connection configuration between the wells and co-construction with existing wells and stations. Similarly, research on star-tree [5] and star-star [6] pipeline networks has become quite mature.A pipeline network also includes facilities such as manifolds, platforms, and compressors. Many studies have investigated oil and gas network models that contain various types of facilities. Rodrigues et al. [7] proposed a 0-1 programming model to determine the connections between FPSO units, manifolds, and facilities as well as the sizes and locations of platforms. In [8, 9], the authors proposed a mixed integer linear programming (MILP) model to optimize the connections between wells and platforms while considering not only the installation time but also the linear drop of the reservoir pressure. Furthermore, Ramos Rosa et al. [10] considered the reservoir dynamics and secondary development to establish an MILP model seeking the maximum net present value (NPV) while simultaneously optimizing the manifold layout, routing, and pipeline diameter. Most of the oil and gas field layout optimization problems mentioned above can be transformed into mixed integer programming (MIP) problems, i.e., programming problems in which some of the independent variables take integer values while others take continuous values. MIP first became an independent branch of mathematical problem-solving when Gomory [11] proposed the cutting-plane method. Generally, MINLP problems can be solved using two types of methods: deterministic algorithms and heuristic algorithms. A deterministic method of solving an MINLP problem is to simplify it to an MILP problem, which can be effectively solved by means of dual simplex or sequential quadratic programming. However, the main difficulty is that appropriate linear approximations must be adopted for nonlinear functions. To this end, Möller [12], Martin et al. [13], Tomasgard [14], and Nørstebø [15] et al. used piecewise approximation to deal with the nonlinear functions in the model. Mikolajková et al. [16] proposed a linearization method for solving MINLP problems. Deterministic methods can guarantee a globally optimal solution when solving an MINLP problem, but they usually require considerable resources to solve large-scale problems.Because of the complexity of MINLP problems, metaheuristic algorithms, including genetic algorithms (GAs) [17] and particle swarm optimization (PSO) [18], are also used to solve problems related to oil and gas pipeline networks. However, the computational burden of heuristic algorithms is usually heavy, making them an inefficient means of solving large-scale problems. In addition, either heuristic methods converge slowly or the quality of their results cannot be guaranteed. Because of the special characteristics of MVs, no related research on optimization problems involving MVs has been reported in recent years; therefore, to address this lack, a corresponding MINLP model and an algorithm for solving it are proposed in this paper.


## 3. Problem Definition and Formulation

### 3.1. Connection Matrix and Assumptions

In general, optimal solutions to various problems can be found through various optimization methods, but the solution found may not necessarily be suitable for the actual conditions and engineering design requirements of a particular situation. For a given oil and gas field, especially an unconventional reservoir with scattered wells, it is necessary to reserve some MV channels for future connections to additional production wells; therefore, either the channel counts of the MVs and the number of MVs need to be fixed in advance or ranges should be defined for these parameters to satisfy the need to reserve valves or channels for future use. If experts define several possible schemes in accordance with the given requirements, the most appropriate solution can be obtained by comparing the results of these schemes, but the quality of a solution chosen in this way cannot be guaranteed; thus, it is necessary to formulate and solve the corresponding optimization problem with constraints on the number of MVs and their channel counts.Suppose that there aren wells that need to be divided into m sets and connected to the corresponding MVs. The MVs will be placed at well sites because no additional land acquisition is allowed; therefore, all well sites are possible alternative locations for MVs, and the number of MVs m is unknown. Let L=lijm×n be the connection matrix, which reflects the connection relationships between the wells and the MVs. If lij=1, then well i is in well set j, which means that well i is connected to MV j. Thus, the remaining elements in the ith column of the matrix are all equal to 0. Here, a matrix Linl1in,l2in,…,lnin,liin∈0,1,i∈1,2,…,n, with dimensions of 1×n is introduced to indicate whether an MV is placed at the ith well site. If liin=1, this means that an MV is placed at the ith well site. If liin=0, this indicates that there is no MV placed at the ith well site. By incorporating Lin into the matrix L, we obtain the matrix L¯, as shown in Figure 2. If liin=1 in the matrix L¯, then the corresponding MV can be identified by finding the row j that contains the element in the ith column of L¯ (below the first row) that has a value of 1; that is, MV j is placed at the ith well site. Thus, ujxjg,yjg and the other wells connected to valve j can be found at the same time. Now, considering the characteristics of MVs, the following assumptions are adopted before solving the mathematical model in this paper:(1)
Generally, the number of channels of an MV ranges from 6 to 16. According to the designer’s requirements, the corresponding MV specification will be given as either a constant or a range before optimization. Thus, the special constraints on the MVs will be specified beforehand.(2)
All wells need to connect directly to an MV because this is the only way in which wells can be selected for metering; there are no pipeline connections between wells.(3)
The possible locations of the MVs are preselected by the engineering company. One of the strategies that is usually adopted is to place the MVs at well sites to reduce land acquisition, which means that all well sites are potential alternative locations for MVs.(4)
In this paper, we use the Euclidean distance to quantify the length of the pipeline between a well and an MV. Additionally, terrain and obstacles are not considered.(5)
The costs of valves and instruments other than the MVs are not considered. These costs are neglected in the calculation since the use of MVs simplifies the metering process.(6)
We assume that the same material and diameter are used for all pipelines. Therefore, factors such as pressure and flow rate are not considered.Figure 2
Matrix of MVs and wells.


### 3.2. MINLP Model

#### 3.2.1. Objective Function

Given the above assumptions, the objective function can be simply expressed as follows:(1)minC=Cpipe+Cv,where Cpipe represents the cost of all pipelines, which is equal to the sum of the costs of the pipelines connected to each MV:(2)Cpipe=∑j=1mCjpipe,where Cjpipe represents the total cost of all pipelines connected to MV j, which can be expressed as follows:(3)Cjpipe=Wp×Lj×Dpipe,j=1,2,…,m,where “×” denotes the matrix dot product operation.The total cost of the pipelines connected to valvej is the dot product of Dpipe, Lj, and Wpw1p,w2p,…,wnp, where Wp is the pipeline cost coefficient matrix, Lj is the jth row of the connection matrix, and Dpipe=d1pipe,d2pipe,…,dnpipe is a matrix consisting of the pipeline lengths from the wells to the valve. oi represents the coordinates of the ith well, i=1,2,…,n. Specifically, oixi,j,yi,j,i=1,2,…,n,j=1,2,…,m, means that the ith well is connected to the jth MV. Then, dipipe represents the pipeline length from well i to the corresponding MV under the assumption that this well is connected to valve j; this pipeline length can be expressed as(4)dipipe=xi,j−xjg2+yi,j−yjg2.u
j
x
j
g
,
y
j
g
∈
O
j represents the coordinates of the jth MV, which are unknown and are determined by solving a subproblem. The detailed solution method will be described in Section 4. Then, the cost of all pipelines can be written as follows:(5)Cpipe=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2.Considering the costs of the MVs themselves, the total cost of all MVs can be expressed as(6)Cv=Wv×B,j=1,2,…,m,where Bβ1,β2,…,βm is a matrix with dimensions of 1×m that is used as a compact representation of the connection matrix, in which βj is a binary variable. If all elements in the jth row of the connection matrix are 0, then βj=0; otherwise, βj=1. Wv is the matrix of the MV cost coefficients. Cv is the dot product of Wv and B. By summing the costs of the pipelines and valves, the objective function is obtained as follows:(7)C=Wp×Lj×Dpipe+Wv×B=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2+∑j=1mwjvβj.


#### 3.2.2. Constraints

The objective function has been established, as shown in equation (7). Now, we will describe the constraints of the model. In this paper, the variable m is constrained by ω¯ and ω¯, as follows:(8)ω¯≤m≤ω¯,where ω¯ is the upper bound on the number of MVs, while ω¯ is the minimum value that m can take; both ω¯ and ω¯ are integers. As the range of m shrinks, the number of MVs will finally be specified; that is, ω¯=m=ω¯. Only one nonzero element is allowed to exist in each column of the connection matrix because each well can be connected to only one MV. This constraint can be expressed as follows:(9)∑j=1mli,j=li,1+li,2⋯li,m=1,i=1,2,…,n,(10)li,j∈0,1.Given that row in the connection matrix represents one MV, the sum of the elements in each row of this matrix cannot be greater than the channel count of the corresponding MV. Since the channel counts of the MVs selected by the designer might differ,Nup is used to represent the maximum MV channel count provided by the manufacturer, and Nr,jup is used to represent the specified channel count of the jth MV:(11)∑i=1nli,j=l1,j+l2,j⋯ln,j≤Njup,j=1,2,…,m.The following constraints should be satisfied:(12)Nr,jup≤Nup,j=1,2,…,m,(13)Njup≤Nup,j=1,2,…,m.IfNr,jup is the channel count of MV j as chosen or required by the designer and ajup is a binary variable that represents whether a specified channel count for the jth valve is actually used to constrain the mathematical model, according to the constraint relaxation degree, then the maximum number of wells that can actually be connected to MV j is defined as(14)Njup=Nup1−ajup+ajupNr,jup,j=1,2,…,m.By combining this equation with equation (11), we obtain(15)∑i=1nli,j≤Nup1−ajup+ajupNr,jup,j=1,2,…,m.Usually, it is not permitted to connect zero wells to a given MV. Therefore, taking 1 as the default lower boundNlb for each row of the matrix, we have(16)∑i=1nli,j=l1,j+l2,j⋯ln,j≥Nlb,j=1,2,…,m.Sometimes, it is unreasonable to connect only a few wells to an MV. In this case, we useNr,jlb to indicate a specified lower bound on the number of wells to which an MV is allowed to be connected. The corresponding constraint can be expressed as(17)∑i=1nli,j=l1,j+l2,j⋯ln,j≥Nr,jlb,j=1,2,…,m.The minimum number of wells that can be connected to MVj is determined by ajlb. ajlb is a binary variable that controls the value of Njlb. If ajlb=1, a specific lower bound is not required and the default is used, whereas if the value of this variable is 0, a specific lower bound is applied. Thus, using the same formulation used for the upper bound, we obtain the following constraint:(18)Njlb=Nlb1−ajlb+ajlbNr,jlb,j=1,2,…,m.By combining equations (16)∼(18), we can write the constraint related to the lower bound as follows:(19)∑i=1nli,j≥Nlb1−ajlb+ajlbNr,jlb,j=1,2,…,m.The following relationships should be satisfied:(20)Nr,jlb≤Njup,j=1,2,…,m,(21)Nr,jlb≤Nr,jup,j=1,2,…,m,(22)Nlb≤Nr,jlb,j=1,2,…,m,(23)Nlb≤Njlb,j=1,2,…,m,(24)Njlb≤Nr,jlb,j=1,2,…,m.Equations (20) and (21) are used to ensure that the lower bound will not exceed the upper bound. Equations (22) and (23) ensure that both the required lower bound and the bound that is actually applied are at least as large as the default lower bound. Finally, equation (24) stipulates that the real lower bound should not exceed the required lower bound. Since the matrix Lin represents the MVs, the sum of the elements of the matrix Lin should satisfy(25)∑i=1nliin=m,liin=0,1.Because only one MV is allowed in each well group, when we subtractLin from any row of the connection matrix L and then add the absolute value of the results together, we will have sumj+m−1 for each row of L. Then, the constraint that only one MV is allowed in each well group can be written as(26)∑i=1nliin−li,j=sumj+m−1,j=1,2,…,m.According to the above equations, the mathematical model can be expressed as(27)C=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2+∑j=1mwjvβj,s.tω¯≤m,m≤ω¯,∑i=1nli,j≤Nup1−ajup+ajupNr,jup,j=1,2,…,m,Nlb1−ajlb+ajlbNr,jlb≤∑i=1nli,j,j=1,2,…,m,Njup≤Nup,j=1,2,…,m,Nr,jup≤Nup,j=1,2,…,m,Nr,jlb≤Njup,j=1,2,…,m,Nr,jlb≤Nr,jup,j=1,2,…,m,Nlb≤Nr,jlb,j=1,2,…,m,Nlb≤Njlb,j=1,2,…,m,Njlb≤Nr,jlb,j=1,2,…,m,∑j=1mli,j=1,i=1,2,…,n,∑i=1nliin=m,∑i=1nliin−li,j=sumj+m−1,j=1,2,…,m,ajup,ajlb,βj,li,j,liin∈0,1.The decision variables used in the mathematical model represent the associations between the wells and the MVs, and the number of valves and their channel counts are treated as constraints. The binary variablesajup and ajlb determine the range of the channel counts of the MVs, while ω¯ and ω¯ determine the range of the number of MVs.


## 3.1. Connection Matrix and Assumptions

In general, optimal solutions to various problems can be found through various optimization methods, but the solution found may not necessarily be suitable for the actual conditions and engineering design requirements of a particular situation. For a given oil and gas field, especially an unconventional reservoir with scattered wells, it is necessary to reserve some MV channels for future connections to additional production wells; therefore, either the channel counts of the MVs and the number of MVs need to be fixed in advance or ranges should be defined for these parameters to satisfy the need to reserve valves or channels for future use. If experts define several possible schemes in accordance with the given requirements, the most appropriate solution can be obtained by comparing the results of these schemes, but the quality of a solution chosen in this way cannot be guaranteed; thus, it is necessary to formulate and solve the corresponding optimization problem with constraints on the number of MVs and their channel counts.Suppose that there aren wells that need to be divided into m sets and connected to the corresponding MVs. The MVs will be placed at well sites because no additional land acquisition is allowed; therefore, all well sites are possible alternative locations for MVs, and the number of MVs m is unknown. Let L=lijm×n be the connection matrix, which reflects the connection relationships between the wells and the MVs. If lij=1, then well i is in well set j, which means that well i is connected to MV j. Thus, the remaining elements in the ith column of the matrix are all equal to 0. Here, a matrix Linl1in,l2in,…,lnin,liin∈0,1,i∈1,2,…,n, with dimensions of 1×n is introduced to indicate whether an MV is placed at the ith well site. If liin=1, this means that an MV is placed at the ith well site. If liin=0, this indicates that there is no MV placed at the ith well site. By incorporating Lin into the matrix L, we obtain the matrix L¯, as shown in Figure 2. If liin=1 in the matrix L¯, then the corresponding MV can be identified by finding the row j that contains the element in the ith column of L¯ (below the first row) that has a value of 1; that is, MV j is placed at the ith well site. Thus, ujxjg,yjg and the other wells connected to valve j can be found at the same time. Now, considering the characteristics of MVs, the following assumptions are adopted before solving the mathematical model in this paper:(1)
Generally, the number of channels of an MV ranges from 6 to 16. According to the designer’s requirements, the corresponding MV specification will be given as either a constant or a range before optimization. Thus, the special constraints on the MVs will be specified beforehand.(2)
All wells need to connect directly to an MV because this is the only way in which wells can be selected for metering; there are no pipeline connections between wells.(3)
The possible locations of the MVs are preselected by the engineering company. One of the strategies that is usually adopted is to place the MVs at well sites to reduce land acquisition, which means that all well sites are potential alternative locations for MVs.(4)
In this paper, we use the Euclidean distance to quantify the length of the pipeline between a well and an MV. Additionally, terrain and obstacles are not considered.(5)
The costs of valves and instruments other than the MVs are not considered. These costs are neglected in the calculation since the use of MVs simplifies the metering process.(6)
We assume that the same material and diameter are used for all pipelines. Therefore, factors such as pressure and flow rate are not considered.Figure 2
Matrix of MVs and wells.


## 3.2. MINLP Model

### 3.2.1. Objective Function

Given the above assumptions, the objective function can be simply expressed as follows:(1)minC=Cpipe+Cv,where Cpipe represents the cost of all pipelines, which is equal to the sum of the costs of the pipelines connected to each MV:(2)Cpipe=∑j=1mCjpipe,where Cjpipe represents the total cost of all pipelines connected to MV j, which can be expressed as follows:(3)Cjpipe=Wp×Lj×Dpipe,j=1,2,…,m,where “×” denotes the matrix dot product operation.The total cost of the pipelines connected to valvej is the dot product of Dpipe, Lj, and Wpw1p,w2p,…,wnp, where Wp is the pipeline cost coefficient matrix, Lj is the jth row of the connection matrix, and Dpipe=d1pipe,d2pipe,…,dnpipe is a matrix consisting of the pipeline lengths from the wells to the valve. oi represents the coordinates of the ith well, i=1,2,…,n. Specifically, oixi,j,yi,j,i=1,2,…,n,j=1,2,…,m, means that the ith well is connected to the jth MV. Then, dipipe represents the pipeline length from well i to the corresponding MV under the assumption that this well is connected to valve j; this pipeline length can be expressed as(4)dipipe=xi,j−xjg2+yi,j−yjg2.u
j
x
j
g
,
y
j
g
∈
O
j represents the coordinates of the jth MV, which are unknown and are determined by solving a subproblem. The detailed solution method will be described in Section 4. Then, the cost of all pipelines can be written as follows:(5)Cpipe=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2.Considering the costs of the MVs themselves, the total cost of all MVs can be expressed as(6)Cv=Wv×B,j=1,2,…,m,where Bβ1,β2,…,βm is a matrix with dimensions of 1×m that is used as a compact representation of the connection matrix, in which βj is a binary variable. If all elements in the jth row of the connection matrix are 0, then βj=0; otherwise, βj=1. Wv is the matrix of the MV cost coefficients. Cv is the dot product of Wv and B. By summing the costs of the pipelines and valves, the objective function is obtained as follows:(7)C=Wp×Lj×Dpipe+Wv×B=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2+∑j=1mwjvβj.


### 3.2.2. Constraints

The objective function has been established, as shown in equation (7). Now, we will describe the constraints of the model. In this paper, the variable m is constrained by ω¯ and ω¯, as follows:(8)ω¯≤m≤ω¯,where ω¯ is the upper bound on the number of MVs, while ω¯ is the minimum value that m can take; both ω¯ and ω¯ are integers. As the range of m shrinks, the number of MVs will finally be specified; that is, ω¯=m=ω¯. Only one nonzero element is allowed to exist in each column of the connection matrix because each well can be connected to only one MV. This constraint can be expressed as follows:(9)∑j=1mli,j=li,1+li,2⋯li,m=1,i=1,2,…,n,(10)li,j∈0,1.Given that row in the connection matrix represents one MV, the sum of the elements in each row of this matrix cannot be greater than the channel count of the corresponding MV. Since the channel counts of the MVs selected by the designer might differ,Nup is used to represent the maximum MV channel count provided by the manufacturer, and Nr,jup is used to represent the specified channel count of the jth MV:(11)∑i=1nli,j=l1,j+l2,j⋯ln,j≤Njup,j=1,2,…,m.The following constraints should be satisfied:(12)Nr,jup≤Nup,j=1,2,…,m,(13)Njup≤Nup,j=1,2,…,m.IfNr,jup is the channel count of MV j as chosen or required by the designer and ajup is a binary variable that represents whether a specified channel count for the jth valve is actually used to constrain the mathematical model, according to the constraint relaxation degree, then the maximum number of wells that can actually be connected to MV j is defined as(14)Njup=Nup1−ajup+ajupNr,jup,j=1,2,…,m.By combining this equation with equation (11), we obtain(15)∑i=1nli,j≤Nup1−ajup+ajupNr,jup,j=1,2,…,m.Usually, it is not permitted to connect zero wells to a given MV. Therefore, taking 1 as the default lower boundNlb for each row of the matrix, we have(16)∑i=1nli,j=l1,j+l2,j⋯ln,j≥Nlb,j=1,2,…,m.Sometimes, it is unreasonable to connect only a few wells to an MV. In this case, we useNr,jlb to indicate a specified lower bound on the number of wells to which an MV is allowed to be connected. The corresponding constraint can be expressed as(17)∑i=1nli,j=l1,j+l2,j⋯ln,j≥Nr,jlb,j=1,2,…,m.The minimum number of wells that can be connected to MVj is determined by ajlb. ajlb is a binary variable that controls the value of Njlb. If ajlb=1, a specific lower bound is not required and the default is used, whereas if the value of this variable is 0, a specific lower bound is applied. Thus, using the same formulation used for the upper bound, we obtain the following constraint:(18)Njlb=Nlb1−ajlb+ajlbNr,jlb,j=1,2,…,m.By combining equations (16)∼(18), we can write the constraint related to the lower bound as follows:(19)∑i=1nli,j≥Nlb1−ajlb+ajlbNr,jlb,j=1,2,…,m.The following relationships should be satisfied:(20)Nr,jlb≤Njup,j=1,2,…,m,(21)Nr,jlb≤Nr,jup,j=1,2,…,m,(22)Nlb≤Nr,jlb,j=1,2,…,m,(23)Nlb≤Njlb,j=1,2,…,m,(24)Njlb≤Nr,jlb,j=1,2,…,m.Equations (20) and (21) are used to ensure that the lower bound will not exceed the upper bound. Equations (22) and (23) ensure that both the required lower bound and the bound that is actually applied are at least as large as the default lower bound. Finally, equation (24) stipulates that the real lower bound should not exceed the required lower bound. Since the matrix Lin represents the MVs, the sum of the elements of the matrix Lin should satisfy(25)∑i=1nliin=m,liin=0,1.Because only one MV is allowed in each well group, when we subtractLin from any row of the connection matrix L and then add the absolute value of the results together, we will have sumj+m−1 for each row of L. Then, the constraint that only one MV is allowed in each well group can be written as(26)∑i=1nliin−li,j=sumj+m−1,j=1,2,…,m.According to the above equations, the mathematical model can be expressed as(27)C=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2+∑j=1mwjvβj,s.tω¯≤m,m≤ω¯,∑i=1nli,j≤Nup1−ajup+ajupNr,jup,j=1,2,…,m,Nlb1−ajlb+ajlbNr,jlb≤∑i=1nli,j,j=1,2,…,m,Njup≤Nup,j=1,2,…,m,Nr,jup≤Nup,j=1,2,…,m,Nr,jlb≤Njup,j=1,2,…,m,Nr,jlb≤Nr,jup,j=1,2,…,m,Nlb≤Nr,jlb,j=1,2,…,m,Nlb≤Njlb,j=1,2,…,m,Njlb≤Nr,jlb,j=1,2,…,m,∑j=1mli,j=1,i=1,2,…,n,∑i=1nliin=m,∑i=1nliin−li,j=sumj+m−1,j=1,2,…,m,ajup,ajlb,βj,li,j,liin∈0,1.The decision variables used in the mathematical model represent the associations between the wells and the MVs, and the number of valves and their channel counts are treated as constraints. The binary variablesajup and ajlb determine the range of the channel counts of the MVs, while ω¯ and ω¯ determine the range of the number of MVs.


## 3.2.1. Objective Function

Given the above assumptions, the objective function can be simply expressed as follows:(1)minC=Cpipe+Cv,where Cpipe represents the cost of all pipelines, which is equal to the sum of the costs of the pipelines connected to each MV:(2)Cpipe=∑j=1mCjpipe,where Cjpipe represents the total cost of all pipelines connected to MV j, which can be expressed as follows:(3)Cjpipe=Wp×Lj×Dpipe,j=1,2,…,m,where “×” denotes the matrix dot product operation.The total cost of the pipelines connected to valvej is the dot product of Dpipe, Lj, and Wpw1p,w2p,…,wnp, where Wp is the pipeline cost coefficient matrix, Lj is the jth row of the connection matrix, and Dpipe=d1pipe,d2pipe,…,dnpipe is a matrix consisting of the pipeline lengths from the wells to the valve. oi represents the coordinates of the ith well, i=1,2,…,n. Specifically, oixi,j,yi,j,i=1,2,…,n,j=1,2,…,m, means that the ith well is connected to the jth MV. Then, dipipe represents the pipeline length from well i to the corresponding MV under the assumption that this well is connected to valve j; this pipeline length can be expressed as(4)dipipe=xi,j−xjg2+yi,j−yjg2.u
j
x
j
g
,
y
j
g
∈
O
j represents the coordinates of the jth MV, which are unknown and are determined by solving a subproblem. The detailed solution method will be described in Section 4. Then, the cost of all pipelines can be written as follows:(5)Cpipe=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2.Considering the costs of the MVs themselves, the total cost of all MVs can be expressed as(6)Cv=Wv×B,j=1,2,…,m,where Bβ1,β2,…,βm is a matrix with dimensions of 1×m that is used as a compact representation of the connection matrix, in which βj is a binary variable. If all elements in the jth row of the connection matrix are 0, then βj=0; otherwise, βj=1. Wv is the matrix of the MV cost coefficients. Cv is the dot product of Wv and B. By summing the costs of the pipelines and valves, the objective function is obtained as follows:(7)C=Wp×Lj×Dpipe+Wv×B=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2+∑j=1mwjvβj.


## 3.2.2. Constraints

The objective function has been established, as shown in equation (7). Now, we will describe the constraints of the model. In this paper, the variable m is constrained by ω¯ and ω¯, as follows:(8)ω¯≤m≤ω¯,where ω¯ is the upper bound on the number of MVs, while ω¯ is the minimum value that m can take; both ω¯ and ω¯ are integers. As the range of m shrinks, the number of MVs will finally be specified; that is, ω¯=m=ω¯. Only one nonzero element is allowed to exist in each column of the connection matrix because each well can be connected to only one MV. This constraint can be expressed as follows:(9)∑j=1mli,j=li,1+li,2⋯li,m=1,i=1,2,…,n,(10)li,j∈0,1.Given that row in the connection matrix represents one MV, the sum of the elements in each row of this matrix cannot be greater than the channel count of the corresponding MV. Since the channel counts of the MVs selected by the designer might differ,Nup is used to represent the maximum MV channel count provided by the manufacturer, and Nr,jup is used to represent the specified channel count of the jth MV:(11)∑i=1nli,j=l1,j+l2,j⋯ln,j≤Njup,j=1,2,…,m.The following constraints should be satisfied:(12)Nr,jup≤Nup,j=1,2,…,m,(13)Njup≤Nup,j=1,2,…,m.IfNr,jup is the channel count of MV j as chosen or required by the designer and ajup is a binary variable that represents whether a specified channel count for the jth valve is actually used to constrain the mathematical model, according to the constraint relaxation degree, then the maximum number of wells that can actually be connected to MV j is defined as(14)Njup=Nup1−ajup+ajupNr,jup,j=1,2,…,m.By combining this equation with equation (11), we obtain(15)∑i=1nli,j≤Nup1−ajup+ajupNr,jup,j=1,2,…,m.Usually, it is not permitted to connect zero wells to a given MV. Therefore, taking 1 as the default lower boundNlb for each row of the matrix, we have(16)∑i=1nli,j=l1,j+l2,j⋯ln,j≥Nlb,j=1,2,…,m.Sometimes, it is unreasonable to connect only a few wells to an MV. In this case, we useNr,jlb to indicate a specified lower bound on the number of wells to which an MV is allowed to be connected. The corresponding constraint can be expressed as(17)∑i=1nli,j=l1,j+l2,j⋯ln,j≥Nr,jlb,j=1,2,…,m.The minimum number of wells that can be connected to MVj is determined by ajlb. ajlb is a binary variable that controls the value of Njlb. If ajlb=1, a specific lower bound is not required and the default is used, whereas if the value of this variable is 0, a specific lower bound is applied. Thus, using the same formulation used for the upper bound, we obtain the following constraint:(18)Njlb=Nlb1−ajlb+ajlbNr,jlb,j=1,2,…,m.By combining equations (16)∼(18), we can write the constraint related to the lower bound as follows:(19)∑i=1nli,j≥Nlb1−ajlb+ajlbNr,jlb,j=1,2,…,m.The following relationships should be satisfied:(20)Nr,jlb≤Njup,j=1,2,…,m,(21)Nr,jlb≤Nr,jup,j=1,2,…,m,(22)Nlb≤Nr,jlb,j=1,2,…,m,(23)Nlb≤Njlb,j=1,2,…,m,(24)Njlb≤Nr,jlb,j=1,2,…,m.Equations (20) and (21) are used to ensure that the lower bound will not exceed the upper bound. Equations (22) and (23) ensure that both the required lower bound and the bound that is actually applied are at least as large as the default lower bound. Finally, equation (24) stipulates that the real lower bound should not exceed the required lower bound. Since the matrix Lin represents the MVs, the sum of the elements of the matrix Lin should satisfy(25)∑i=1nliin=m,liin=0,1.Because only one MV is allowed in each well group, when we subtractLin from any row of the connection matrix L and then add the absolute value of the results together, we will have sumj+m−1 for each row of L. Then, the constraint that only one MV is allowed in each well group can be written as(26)∑i=1nliin−li,j=sumj+m−1,j=1,2,…,m.According to the above equations, the mathematical model can be expressed as(27)C=∑j=1m∑i=1nwipli,jxi,j−xjg2+yi,j−yjg2+∑j=1mwjvβj,s.tω¯≤m,m≤ω¯,∑i=1nli,j≤Nup1−ajup+ajupNr,jup,j=1,2,…,m,Nlb1−ajlb+ajlbNr,jlb≤∑i=1nli,j,j=1,2,…,m,Njup≤Nup,j=1,2,…,m,Nr,jup≤Nup,j=1,2,…,m,Nr,jlb≤Njup,j=1,2,…,m,Nr,jlb≤Nr,jup,j=1,2,…,m,Nlb≤Nr,jlb,j=1,2,…,m,Nlb≤Njlb,j=1,2,…,m,Njlb≤Nr,jlb,j=1,2,…,m,∑j=1mli,j=1,i=1,2,…,n,∑i=1nliin=m,∑i=1nliin−li,j=sumj+m−1,j=1,2,…,m,ajup,ajlb,βj,li,j,liin∈0,1.The decision variables used in the mathematical model represent the associations between the wells and the MVs, and the number of valves and their channel counts are treated as constraints. The binary variablesajup and ajlb determine the range of the channel counts of the MVs, while ω¯ and ω¯ determine the range of the number of MVs.


## 4. The Proposed Algorithm for Solving the MINLP Problem

The underlying principle of methods such as steepest gradient descent [19] is to find the direction in which the objective function decreases the fastest at the present point, take a step in that direction, and then repeat this process, thus iteratively reducing the value of the objective function. We use the same idea to search for the optimal solution in a discrete space. In this paper, the mathematical model is described in the matrix form; thus, the step size and the number of steps must be expressed in terms of a measure of the distance between two matrices. If we treat the whole matrix as a point, then we can move one nonzero element in the matrix to obtain a new matrix, and this new matrix can be treated as a point adjacent to the original point; that is, the new matrix is one step away from the original matrix. All the points obtained by changing one individual element are the points closest to the current point. As an example, let us consider a connection matrix that is equivalent to the 5 × 5 unit matrix; then, we can obtain a foresight matrix corresponding to this connection matrix, as shown in Figure 3. The square columns represent the objective function values, where each row of the matrix represents a well and each column represents an MV. Then, the values in the jth column of the matrix represent the values of the objective function after the corresponding nonzero elements are moved to the jth valve; that is, the square column in the ith row and jth column represents the objective function value after well i is connected to MV j.Figure 3
Foresight matrix of the unit matrix.If the initial connection matrix is a unit matrix, then initially, the values on the diagonal of the foresight matrix are equal to the objective function value of the original matrix because the values on the diagonal of the foresight matrix, represented by square columns with stripes, are obtained only when no element is moved. For a case in which the constraints are not satisfied after a single nonzero element is moved to another position, the corresponding objective function value is represented by a dark gray square column; all such objective function values are replaced by a sufficiently large penalty valueM to ensure that the algorithm will not select such a point. White columns are used to represent the feasible domain, meaning that these positions in the matrix represent points that satisfy the constraints of the mathematical model. The minimum value in the feasible domain, denoted by vn+1, is found as the result of the current iteration, as shown in Figure 3. As shown in the figure, the minimum objective function value is found in the 4th row and 1st column; therefore, we set the element in the 4th row and 1st column to 1 and set the remaining elements in the 4th row to 0. In this way, we can find the path that causes the objective function value to drop the fastest from the current point, thus completing one single-step iteration. When all values in the feasible domain are greater than or equal to the current value, the current point is considered to be a locally optimal solution.We can always construct a foresight matrix such as that shown in Figure3 after every step and then find the point that minimizes the objective function for the next step until no value less than the current value of the objective function is found. To obtain a better solution once the current point has fallen into a local optimum, we can consider increasing the step size, thus enabling the algorithm to escape from the local optimum and find a better feasible solution after further iterations.


### 4.1. Solving the Subproblem

All wells connected to MVj are assigned to the same well set. The MV will be placed at one of the wells in this set such that the total weighted length from the wells to the MV is as short as possible, as expressed in equation (23). Thus, a subproblem must be solved to find the placement of valve j that results in the shortest pipeline length in well group j:(28)minSj∑i=1nli,jdipipe,j=1,2,…,m.This subproblem can be solved using a simple algorithm. The coordinates of the wells connected to valvej are represented by xi,j,yi,j,i=1,…,numj,j=1,…,m. Because the location of MV j is unknown, ujxjg,yjg∈Oj is assumed to represent the valve coordinates that minimize the sum of the pipeline lengths connected to valve j. The valve will be placed at each well site in turn; then, we can compare the results to find the optimal placement position for the valve.


### 4.2. Details of the Algorithm

The algorithm for solving the MINLP problem obtains the solution through iteration. The iterative process mainly comprises 4 steps; this process is depicted in Figure4 and described in detail below.Figure 4
Algorithm flow diagram.


#### 4.2.1. Step 1: Initialization

If the constraint on the number of MVsm is relaxed, m can range from 1 to n. In the extreme case, a unit matrix is taken as the initial feasible solution, meaning that each MV is initially assumed to be connected to only one well. Generally, if 1≤ω¯≤m≤ω¯≤n, we can initially set m=ω¯, and the wells will be sequentially assigned to the n MVs in accordance with the constraint conditions.


#### 4.2.2. Step 2: Iteration

The initial costs can be obtained by substituting the initial distribution matrix into the objective function. Then, a blank matrix with dimensions ofm×n is initialized as the foresight matrix.Starting from the first column, each nonzero element 1 in the connection matrix is moved up or down by one row; each movement will change the distribution, and the objective function value will change as the connection matrix changes. The corresponding objective function value is calculated for every movement, and the new objective function value is filled into the corresponding position of the foresight matrix.In the example shown in Figure5, element 1 in the first row and first column is moved to the second row and first column to obtain a new connection relationship. After substituting this new connection relationship into the objective function, the value v1,2 corresponding to the new connection relationship is obtained. v1,2 is then filled into the second row and first column of the foresight matrix, corresponding to the new position of the nonzero element in the connection matrix. If the constraints are not satisfied by this new connection relationship, a sufficiently large penalty value M is placed in the corresponding position instead. Once the foresight matrix has been completely filled, the path of fastest descent can be determined. This process is illustrated in Figure 5.Figure 5
Process of filling in the foresight matrix.According to the constraints, the nonzero elements in the connection matrix can move only up and down; they cannot move left or right. Since each column of the matrix represents a well, when a nonzero element is moved from thejth row to the kth row in column i, this represents that this well is disconnected from MV j and connected to MV k, and the value of the element in the ith row and jth column returns to 0.In this way, each 1 value in the connection matrix is moved from the top row to the bottom. Each time a 1 value is moved to a new position, a new objective function value is generated, filling in the corresponding position in the foresight matrix. When the new matrix represents a connection relationship that does not satisfy the constraints, a sufficiently large penalty valueM is filled into the foresight matrix instead. Once the m×n foresight matrix has been completed, we then search for the minimum value in this matrix that is smaller than the current value. For example, if vi,j is the smallest value in the foresight matrix, then for the next iteration, 1 in the ith column of the connection matrix is moved to the jth row, and the value at the original position is set to 0; thus, one single-step iteration is completed. By repeating this process, the current value of the objective function can be iteratively reduced, and new connection relationships can be formed. Figure 5 illustrates this iterative process for a 5 × 5 foresight matrix.If all nonzero elements in a given row of the connection matrix leave that row by moving up or down, this means that the corresponding MV does not exist; in this case, the empty row should be removed from the connection matrix, decreasing the feasible domain. In this way, the number of MVs can be found. Specifically, if the elements of thejth row are all 0, then the jth row should be removed from the connection matrix, as shown in Figure 6.Figure 6
Removal of a row that does not contain any wells.


#### 4.2.3. Step 3: Increasing the Step Size

If a lower cost cannot be obtained by moving a single element in the connection matrix up or down, this means that the solution has become stuck in a local optimum. Therefore, the step size must be gradually increased to allow the algorithm to escape from this local optimum. For example, suppose that the matrixL is initially set as the unit matrix and that the initial step size is 1. When the objective function value cannot continue to decrease after a certain number of iterations, then increasing the step size can help to find a new connection relationship that will allow the objective function to continue to decrease. More specifically, suppose that the ith row of the matrix L contains two nonzero elements, meaning that MV i is connected to two wells, while the jth row contains three nonzero elements. In this case, both nonzero elements in the ith row are added to the jth row, setting all elements in the ith row to zero, and the new total cost vi,j under the new connection relationship is obtained. This cost value is filled into the ith row and jth column of the foresight matrix. Then, the matrix is restored to its previous state, and the algorithm proceeds to fill in the other elements in the foresight matrix in the same way. Each element in the foresight matrix is an objective function value, and these values are used as the basis for iteration.In the foresight matrix for multistep iteration, the values on the diagonal should be equivalent to the objective function value when no wells are moved because for each element on the diagonal, the row and column correspond to the same well; therefore, the penalty valueM is filled into the positions on the diagonal of the matrix. Moreover, the matrix elements below the diagonal to the left are symmetric to the matrix elements above the diagonal to the right; therefore, all matrix elements in either the lower left or the upper right can be directly set equal to M to avoid redundant calculations. As before, if the constraints are not satisfied after the elements of the connection matrix are moved, the penalty value M is also directly filled into the corresponding position. After the foresight matrix has been fully filled in, the smallest value in this matrix can be found. If this minimum value is smaller than the cost obtained in the previous iteration, then we can find the row and column corresponding to this minimum value, add the elements of the ith row to the jth row, and remove the ith row, whose elements are now all 0. At this time, the number of MVs m is set to m minus 1.If the remaining channels of an MV do not allow connection to any further wells, the wells will be sequentially assigned to other MVs. This iterative process will repeat until the objective function value can no longer be reduced by such multistep iterations. Then, the algorithm will either return to Step 1 or proceed to Step 4.


#### 4.2.4. Step 4: Adjustment

If there is no MV whose channels are fully occupied or for which the number of wells connected to that valve has reached the lower boundNjlb, then the optimal solution can be obtained through single-step or multistep iteration. However, according to the constraints, if an MV is fully occupied, then it cannot be connected to any new wells, whereas if a valve has reached its lower bound Njlb, then no more wells can be moved out of the corresponding row. In this case, the iterative process can be continued as follows:(1)
The step size can be increased to 2 to test whether the objective function value will continue to decrease. First, we look for the fully occupied MV. Suppose thatnumj wells are connected to this MV; then, the remaining n−numj wells will each be moved to this MV in turn, replacing one of the existing wells, which is moved to another MV that meets the constraints. By repeating this process, the new values of the objective function are obtained and written into the foresight matrix based on the corresponding well numbers i and j. Finally, the path of steepest descent can be determined in the same way as before.(2)
For an MVj that has been fully connected, we look for an MV t with more channels than MV j, where bj denotes the number of channels of MV j. Then, bj<bt, and the total number of wells connected to each valve satisfies sumt≤sumj. Since sumj=bj, bj≥sumt must be satisfied to allow the jth and tth valves to be exchanged; that is, the following conditions should be satisfied simultaneously:(29)bj<bt,bj≥sumt.If these conditions are satisfied, then the channel counts of thejth and tth MVs are exchanged, and we continue iterating steps 2 and 3 mentioned above. However, if i and t do not satisfy bisum<bt, bi<bt, and sumb,t≤sumb,t, then the penalty value M is written into the corresponding position of the foresight matrix. The process above is repeated; if, after the Nth iteration, the objective function value cannot continue to decrease, meaning that vN+1min−vNmin=0, then the current objective function value corresponds to the optimal solution.


## 4.1. Solving the Subproblem

All wells connected to MVj are assigned to the same well set. The MV will be placed at one of the wells in this set such that the total weighted length from the wells to the MV is as short as possible, as expressed in equation (23). Thus, a subproblem must be solved to find the placement of valve j that results in the shortest pipeline length in well group j:(28)minSj∑i=1nli,jdipipe,j=1,2,…,m.This subproblem can be solved using a simple algorithm. The coordinates of the wells connected to valvej are represented by xi,j,yi,j,i=1,…,numj,j=1,…,m. Because the location of MV j is unknown, ujxjg,yjg∈Oj is assumed to represent the valve coordinates that minimize the sum of the pipeline lengths connected to valve j. The valve will be placed at each well site in turn; then, we can compare the results to find the optimal placement position for the valve.


## 4.2. Details of the Algorithm

The algorithm for solving the MINLP problem obtains the solution through iteration. The iterative process mainly comprises 4 steps; this process is depicted in Figure4 and described in detail below.Figure 4
Algorithm flow diagram.


### 4.2.1. Step 1: Initialization

If the constraint on the number of MVsm is relaxed, m can range from 1 to n. In the extreme case, a unit matrix is taken as the initial feasible solution, meaning that each MV is initially assumed to be connected to only one well. Generally, if 1≤ω¯≤m≤ω¯≤n, we can initially set m=ω¯, and the wells will be sequentially assigned to the n MVs in accordance with the constraint conditions.


### 4.2.2. Step 2: Iteration

The initial costs can be obtained by substituting the initial distribution matrix into the objective function. Then, a blank matrix with dimensions ofm×n is initialized as the foresight matrix.Starting from the first column, each nonzero element 1 in the connection matrix is moved up or down by one row; each movement will change the distribution, and the objective function value will change as the connection matrix changes. The corresponding objective function value is calculated for every movement, and the new objective function value is filled into the corresponding position of the foresight matrix.In the example shown in Figure5, element 1 in the first row and first column is moved to the second row and first column to obtain a new connection relationship. After substituting this new connection relationship into the objective function, the value v1,2 corresponding to the new connection relationship is obtained. v1,2 is then filled into the second row and first column of the foresight matrix, corresponding to the new position of the nonzero element in the connection matrix. If the constraints are not satisfied by this new connection relationship, a sufficiently large penalty value M is placed in the corresponding position instead. Once the foresight matrix has been completely filled, the path of fastest descent can be determined. This process is illustrated in Figure 5.Figure 5
Process of filling in the foresight matrix.According to the constraints, the nonzero elements in the connection matrix can move only up and down; they cannot move left or right. Since each column of the matrix represents a well, when a nonzero element is moved from thejth row to the kth row in column i, this represents that this well is disconnected from MV j and connected to MV k, and the value of the element in the ith row and jth column returns to 0.In this way, each 1 value in the connection matrix is moved from the top row to the bottom. Each time a 1 value is moved to a new position, a new objective function value is generated, filling in the corresponding position in the foresight matrix. When the new matrix represents a connection relationship that does not satisfy the constraints, a sufficiently large penalty valueM is filled into the foresight matrix instead. Once the m×n foresight matrix has been completed, we then search for the minimum value in this matrix that is smaller than the current value. For example, if vi,j is the smallest value in the foresight matrix, then for the next iteration, 1 in the ith column of the connection matrix is moved to the jth row, and the value at the original position is set to 0; thus, one single-step iteration is completed. By repeating this process, the current value of the objective function can be iteratively reduced, and new connection relationships can be formed. Figure 5 illustrates this iterative process for a 5 × 5 foresight matrix.If all nonzero elements in a given row of the connection matrix leave that row by moving up or down, this means that the corresponding MV does not exist; in this case, the empty row should be removed from the connection matrix, decreasing the feasible domain. In this way, the number of MVs can be found. Specifically, if the elements of thejth row are all 0, then the jth row should be removed from the connection matrix, as shown in Figure 6.Figure 6
Removal of a row that does not contain any wells.


### 4.2.3. Step 3: Increasing the Step Size

If a lower cost cannot be obtained by moving a single element in the connection matrix up or down, this means that the solution has become stuck in a local optimum. Therefore, the step size must be gradually increased to allow the algorithm to escape from this local optimum. For example, suppose that the matrixL is initially set as the unit matrix and that the initial step size is 1. When the objective function value cannot continue to decrease after a certain number of iterations, then increasing the step size can help to find a new connection relationship that will allow the objective function to continue to decrease. More specifically, suppose that the ith row of the matrix L contains two nonzero elements, meaning that MV i is connected to two wells, while the jth row contains three nonzero elements. In this case, both nonzero elements in the ith row are added to the jth row, setting all elements in the ith row to zero, and the new total cost vi,j under the new connection relationship is obtained. This cost value is filled into the ith row and jth column of the foresight matrix. Then, the matrix is restored to its previous state, and the algorithm proceeds to fill in the other elements in the foresight matrix in the same way. Each element in the foresight matrix is an objective function value, and these values are used as the basis for iteration.In the foresight matrix for multistep iteration, the values on the diagonal should be equivalent to the objective function value when no wells are moved because for each element on the diagonal, the row and column correspond to the same well; therefore, the penalty valueM is filled into the positions on the diagonal of the matrix. Moreover, the matrix elements below the diagonal to the left are symmetric to the matrix elements above the diagonal to the right; therefore, all matrix elements in either the lower left or the upper right can be directly set equal to M to avoid redundant calculations. As before, if the constraints are not satisfied after the elements of the connection matrix are moved, the penalty value M is also directly filled into the corresponding position. After the foresight matrix has been fully filled in, the smallest value in this matrix can be found. If this minimum value is smaller than the cost obtained in the previous iteration, then we can find the row and column corresponding to this minimum value, add the elements of the ith row to the jth row, and remove the ith row, whose elements are now all 0. At this time, the number of MVs m is set to m minus 1.If the remaining channels of an MV do not allow connection to any further wells, the wells will be sequentially assigned to other MVs. This iterative process will repeat until the objective function value can no longer be reduced by such multistep iterations. Then, the algorithm will either return to Step 1 or proceed to Step 4.


### 4.2.4. Step 4: Adjustment

If there is no MV whose channels are fully occupied or for which the number of wells connected to that valve has reached the lower boundNjlb, then the optimal solution can be obtained through single-step or multistep iteration. However, according to the constraints, if an MV is fully occupied, then it cannot be connected to any new wells, whereas if a valve has reached its lower bound Njlb, then no more wells can be moved out of the corresponding row. In this case, the iterative process can be continued as follows:(1)
The step size can be increased to 2 to test whether the objective function value will continue to decrease. First, we look for the fully occupied MV. Suppose thatnumj wells are connected to this MV; then, the remaining n−numj wells will each be moved to this MV in turn, replacing one of the existing wells, which is moved to another MV that meets the constraints. By repeating this process, the new values of the objective function are obtained and written into the foresight matrix based on the corresponding well numbers i and j. Finally, the path of steepest descent can be determined in the same way as before.(2)
For an MVj that has been fully connected, we look for an MV t with more channels than MV j, where bj denotes the number of channels of MV j. Then, bj<bt, and the total number of wells connected to each valve satisfies sumt≤sumj. Since sumj=bj, bj≥sumt must be satisfied to allow the jth and tth valves to be exchanged; that is, the following conditions should be satisfied simultaneously:(29)bj<bt,bj≥sumt.If these conditions are satisfied, then the channel counts of thejth and tth MVs are exchanged, and we continue iterating steps 2 and 3 mentioned above. However, if i and t do not satisfy bisum<bt, bi<bt, and sumb,t≤sumb,t, then the penalty value M is written into the corresponding position of the foresight matrix. The process above is repeated; if, after the Nth iteration, the objective function value cannot continue to decrease, meaning that vN+1min−vNmin=0, then the current objective function value corresponds to the optimal solution.


## 4.2.1. Step 1: Initialization

If the constraint on the number of MVsm is relaxed, m can range from 1 to n. In the extreme case, a unit matrix is taken as the initial feasible solution, meaning that each MV is initially assumed to be connected to only one well. Generally, if 1≤ω¯≤m≤ω¯≤n, we can initially set m=ω¯, and the wells will be sequentially assigned to the n MVs in accordance with the constraint conditions.


## 4.2.2. Step 2: Iteration

The initial costs can be obtained by substituting the initial distribution matrix into the objective function. Then, a blank matrix with dimensions ofm×n is initialized as the foresight matrix.Starting from the first column, each nonzero element 1 in the connection matrix is moved up or down by one row; each movement will change the distribution, and the objective function value will change as the connection matrix changes. The corresponding objective function value is calculated for every movement, and the new objective function value is filled into the corresponding position of the foresight matrix.In the example shown in Figure5, element 1 in the first row and first column is moved to the second row and first column to obtain a new connection relationship. After substituting this new connection relationship into the objective function, the value v1,2 corresponding to the new connection relationship is obtained. v1,2 is then filled into the second row and first column of the foresight matrix, corresponding to the new position of the nonzero element in the connection matrix. If the constraints are not satisfied by this new connection relationship, a sufficiently large penalty value M is placed in the corresponding position instead. Once the foresight matrix has been completely filled, the path of fastest descent can be determined. This process is illustrated in Figure 5.Figure 5
Process of filling in the foresight matrix.According to the constraints, the nonzero elements in the connection matrix can move only up and down; they cannot move left or right. Since each column of the matrix represents a well, when a nonzero element is moved from thejth row to the kth row in column i, this represents that this well is disconnected from MV j and connected to MV k, and the value of the element in the ith row and jth column returns to 0.In this way, each 1 value in the connection matrix is moved from the top row to the bottom. Each time a 1 value is moved to a new position, a new objective function value is generated, filling in the corresponding position in the foresight matrix. When the new matrix represents a connection relationship that does not satisfy the constraints, a sufficiently large penalty valueM is filled into the foresight matrix instead. Once the m×n foresight matrix has been completed, we then search for the minimum value in this matrix that is smaller than the current value. For example, if vi,j is the smallest value in the foresight matrix, then for the next iteration, 1 in the ith column of the connection matrix is moved to the jth row, and the value at the original position is set to 0; thus, one single-step iteration is completed. By repeating this process, the current value of the objective function can be iteratively reduced, and new connection relationships can be formed. Figure 5 illustrates this iterative process for a 5 × 5 foresight matrix.If all nonzero elements in a given row of the connection matrix leave that row by moving up or down, this means that the corresponding MV does not exist; in this case, the empty row should be removed from the connection matrix, decreasing the feasible domain. In this way, the number of MVs can be found. Specifically, if the elements of thejth row are all 0, then the jth row should be removed from the connection matrix, as shown in Figure 6.Figure 6
Removal of a row that does not contain any wells.


## 4.2.3. Step 3: Increasing the Step Size

If a lower cost cannot be obtained by moving a single element in the connection matrix up or down, this means that the solution has become stuck in a local optimum. Therefore, the step size must be gradually increased to allow the algorithm to escape from this local optimum. For example, suppose that the matrixL is initially set as the unit matrix and that the initial step size is 1. When the objective function value cannot continue to decrease after a certain number of iterations, then increasing the step size can help to find a new connection relationship that will allow the objective function to continue to decrease. More specifically, suppose that the ith row of the matrix L contains two nonzero elements, meaning that MV i is connected to two wells, while the jth row contains three nonzero elements. In this case, both nonzero elements in the ith row are added to the jth row, setting all elements in the ith row to zero, and the new total cost vi,j under the new connection relationship is obtained. This cost value is filled into the ith row and jth column of the foresight matrix. Then, the matrix is restored to its previous state, and the algorithm proceeds to fill in the other elements in the foresight matrix in the same way. Each element in the foresight matrix is an objective function value, and these values are used as the basis for iteration.In the foresight matrix for multistep iteration, the values on the diagonal should be equivalent to the objective function value when no wells are moved because for each element on the diagonal, the row and column correspond to the same well; therefore, the penalty valueM is filled into the positions on the diagonal of the matrix. Moreover, the matrix elements below the diagonal to the left are symmetric to the matrix elements above the diagonal to the right; therefore, all matrix elements in either the lower left or the upper right can be directly set equal to M to avoid redundant calculations. As before, if the constraints are not satisfied after the elements of the connection matrix are moved, the penalty value M is also directly filled into the corresponding position. After the foresight matrix has been fully filled in, the smallest value in this matrix can be found. If this minimum value is smaller than the cost obtained in the previous iteration, then we can find the row and column corresponding to this minimum value, add the elements of the ith row to the jth row, and remove the ith row, whose elements are now all 0. At this time, the number of MVs m is set to m minus 1.If the remaining channels of an MV do not allow connection to any further wells, the wells will be sequentially assigned to other MVs. This iterative process will repeat until the objective function value can no longer be reduced by such multistep iterations. Then, the algorithm will either return to Step 1 or proceed to Step 4.


## 4.2.4. Step 4: Adjustment

If there is no MV whose channels are fully occupied or for which the number of wells connected to that valve has reached the lower boundNjlb, then the optimal solution can be obtained through single-step or multistep iteration. However, according to the constraints, if an MV is fully occupied, then it cannot be connected to any new wells, whereas if a valve has reached its lower bound Njlb, then no more wells can be moved out of the corresponding row. In this case, the iterative process can be continued as follows:(1)
The step size can be increased to 2 to test whether the objective function value will continue to decrease. First, we look for the fully occupied MV. Suppose thatnumj wells are connected to this MV; then, the remaining n−numj wells will each be moved to this MV in turn, replacing one of the existing wells, which is moved to another MV that meets the constraints. By repeating this process, the new values of the objective function are obtained and written into the foresight matrix based on the corresponding well numbers i and j. Finally, the path of steepest descent can be determined in the same way as before.(2)
For an MVj that has been fully connected, we look for an MV t with more channels than MV j, where bj denotes the number of channels of MV j. Then, bj<bt, and the total number of wells connected to each valve satisfies sumt≤sumj. Since sumj=bj, bj≥sumt must be satisfied to allow the jth and tth valves to be exchanged; that is, the following conditions should be satisfied simultaneously:(29)bj<bt,bj≥sumt.If these conditions are satisfied, then the channel counts of thejth and tth MVs are exchanged, and we continue iterating steps 2 and 3 mentioned above. However, if i and t do not satisfy bisum<bt, bi<bt, and sumb,t≤sumb,t, then the penalty value M is written into the corresponding position of the foresight matrix. The process above is repeated; if, after the Nth iteration, the objective function value cannot continue to decrease, meaning that vN+1min−vNmin=0, then the current objective function value corresponds to the optimal solution.


## 5. Analysis of Examples

### 5.1. Verification of the Validity of the Algorithm

First, a simple example is presented to demonstrate that the proposed algorithm can effectively solve the problem of interest. The associations between the wells and the MVs and the number of MVs are taken as variables, and their values are found by the algorithm via iteration; in particular, during the iterative process, multistep iteration can help the objective function value to escape from a local optimum. Let us consider a simple example for which we can make a reasonable prediction of the expected solution. Then, we can apply the proposed algorithm to solve the example and compare its solution with our expected result to prove that the algorithm can effectively solve the presented MINLP model.For all examples in this paper, we consider steel pipes with an inner diameter of 53 mm and a wall thickness of 3.5 mm as the pipes linking the wells and valves; thus, the cost coefficientwip is a constant. Moreover, the cost of an MV wjv is a known constant once the vender is decided. In this paper, the cost of an MV is treated as an average value that does not vary with the number of channels. Here, we consider an example of a field with 26 wells (Table 1; W10∼W35), for which the channel counts of the MVs and the number of MVs are not limited; therefore, the number of channels of each MV is uniformly set to 16, which is the maximum number of channels per MV that can be provided by general manufacturers.We first perform single-step iteration, during which the number of MVs will usually stay at a relatively large value, as shown in Figure7. When the identity matrix is used as the initial matrix, the number of MVs cannot be decreased below 8 via this process. However, increasing the step size can lead to further reduction in the value of the objective function, demonstrating that the reason that the objective function cannot be further reduced via single-step iteration is that the step size of each advance is too small. Thus, the objective function becomes stuck in a local optimum, and it is necessary to increase the step size to escape from that local optimum to allow the objective function to continue to decrease. The objective function decreases almost linearly during the initial 18 steps because the main factor affecting the decrease in the objective function value during the initial iterations is the number of MVs rather than the pipeline cost. However, the value of the objective function no longer decreases once it drops into a local optimum. After the 18th iteration, the step size must be increased to continue to reduce the objective function value. The single-step iteration process converges prematurely after the 18th iteration; thus, the findings demonstrate that the strategy of increasing the step size can efficiently help the objective function to escape from a local optimum.Figure 7
Comparison of two different step sizes.The result of optimization with an increase in step size is shown in Figure8. The hollow square symbols represent the positions of the MVs, the triangular symbols represent the wells connected to MV 1, and the circular symbols represent the wells connected to MV 2. For the 26 wells considered in this calculation, the algorithm completely divides them into two groups, and the solution satisfies all the constraints, which is obviously consistent with our expectations. This example proves that the proposed algorithm can be applied to solve the problem of interest and obtain reasonable optimization results.Figure 8
Results of optimization with an increase in step size.


### 5.2. Convergence of the Algorithm

Usually, we can judge whether an algorithm is valid according to whether the value of the objective function converges during the iterative process. If the algorithm is divergent, no matter how many times it is iterated, the value of the objective function will always fluctuate and will not approach a constant. In contrast, convergence of the algorithm means that the objective function value will tend toward stability after some numbers of iterations. In summary, convergence is a necessary condition for the algorithm to be valid.Next, various examples are presented to analyze the convergence of the algorithm. First, we consider a field with 35 wells (Table2; W1∼W35) as an example. MVs with four different channel counts are specified: N1up=14,N2up=12,N3up=10,andN4up=8. Thus, the number of channels of each MV is fixed in advance, the limitations on the numbers of pipelines that are allowed to be connected to different valves are different, and the number of MVs is also given, i.e., m = 4.Figure9 shows the iterative process for these 35 wells, which includes single-step iteration, double-step iteration, and redistribution of different kinds of valves. The iterative process still converges, and the optimal solution is obtained after the 21st iteration. After 14 iterations, the redistribution of different valve types is necessary to allow the objective function value to continue to decrease. Finally, a difference of 0 between two consecutive steps indicates that the final solution has been obtained.Figure 9
Iterative process for 35 wells under a fixed number of MVs with fixed channel counts.As seen from Figure9, the objective function value gradually tends toward stability and converges as iteration continues. The result is shown in Figure 10. Red circular symbols, blue triangular symbols, purple diamond symbols, and yellow triangular symbols are used to represent the individual wells connected to MVs 1, 2, 3, and 4, respectively, and square symbols are used to represent the MVs. The initial channel counts of the MVs corresponding to these different well groups are 14, 12, 10, and 8, respectively, while the final channel counts are 12, 10, 14, and 8, in sequence. In the case of tight constraints, the minimum value of the objective function can be obtained by using the proposed optimization algorithm.Figure 10
Optimization results for 35 wells.In addition to the example above, an example with 40 wells (Table1; W69∼108) and various channel counts is considered to prove the convergence of the algorithm; the results are shown in Figures 11(a) and 11(b). It can be seen from the figures that, during the initial iterations, the objective function value decreases at a relatively fast rate. Then, the rate of decrease slowly declines as iteration proceeds. The algorithm eventually converges when vN−1min−vNmin=0.Figure 11
Convergence examples under different constraints: (a) 6 valves, 40∼48 channels; (b) 8 valves, 40∼48 channels.
(a)
(b)Again, the value of the objective function gradually converges and tends toward stability. In Figure11(a), the iterative processes with a total of 46 and 48 channels are identical, indicating that once the number of channels exceeds a certain threshold, a further increase in the number of channels has no effect on the iterative process or the result. Additionally, although the iterative process for an initial feasible solution with 40 channels is different from the other 4 cases, the final result is not significantly different from the others, indicating that the initial solution does not have a great influence on the result.The iterative processes for total channel counts of 46 and 48 are also fundamentally similar. Once the total number of channels exceeds 46, the iterative processes are basically the same. The results show that the algorithm presented above can stably converge. Therefore, these examples demonstrate that even under different constraints and in different scenarios, the algorithm can drive the objective function value toward stable convergence.


### 5.3. Parameter Analysis

In general, the designer of a pipeline system should consider the possibility of reserving channels for new wells in the future along with other factors when specifying the channel counts of the MVs and the number of MVs. However, different combinations of MVs with different numbers of channels might affect the results; therefore, we will analyze the effects of different channel counts and different numbers of MVs on the results of our algorithm. We consider a field with 40 wells (Table1; W69∼108) as our example here, varying the constraints on the numbers of channels and valves to obtain different results. The results are shown in Table 6, from which we can see the influence of these different constraint parameters on the cost of the pipeline system.Table 6
Costs with various numbers of channels and MVs.
Total channels
Cost (¥)
Number of MVs
Individual channel counts
40
1.424 × 106
3
16, 16, 8
42
1.396 × 106
3
16, 16, 10
44
1.396 × 106
3
16, 14, 14
46
1.396 × 106
3
16, 16, 14
48
1.396 × 106
3
16, 16, 16
50
1.304 × 106
3
16, 16, 18
40
1.306 × 106
4
10, 10, 10, 10
42
1.287 × 106
4
12, 10, 10, 10
44
1.287 × 106
4
12, 12, 10, 10
46
1.287 × 106
4
12, 12, 12, 10
48
1.287 × 106
4
12, 12, 12, 12
50
1.287 × 106
4
14, 12, 12, 12
40
1.647 × 106
5
8, 8, 8, 8, 8
42
1.520 × 106
5
8, 8, 8, 8, 10
44
1.426 × 106
5
8, 8, 8, 10, 10
46
1.426 × 106
5
8, 8, 10, 10, 10
48
1.426 × 106
5
8, 8, 10, 10, 12
50
1.426 × 106
5
8, 8, 10, 12, 12
40
1.641 × 106
6
7, 6, 7, 7, 7, 6
42
1.619 × 106
6
8, 8, 6, 6, 6, 8
44
1.615 × 106
6
8, 8, 8, 8, 6, 6
46
1.592 × 106
6
8, 8, 8, 10, 6, 6
48
1.592 × 106
6
10, 8, 8, 10, 6, 6
50
1.592 × 106
6
10, 8, 8, 8, 8, 8
40
1.893 × 106
7
6, 6, 6, 6, 6, 4, 6
42
1.870 × 106
7
6, 6, 6, 6, 6, 6, 6
44
1.775 × 106
7
6, 8, 6, 6, 6, 6, 6
46
1.766 × 106
7
8, 8, 6, 6, 6, 6, 6
48
1.766 × 106
7
6, 8, 6, 8, 8, 6, 6
50
1.766 × 106
7
8, 8, 6, 8, 6, 6, 6
40
1.957 × 106
8
5, 5, 5, 5, 5, 5, 5, 5
42
1.950 × 106
8
6, 6, 5, 5, 5, 5, 5, 5
44
1.941 × 106
8
6, 6, 5, 5, 6, 6, 5, 5
46
1.941 × 106
8
6, 6, 5, 5, 6, 6, 6, 6
48
1.941 × 106
8
6, 6, 5, 5, 6, 6, 8, 6
50
1.919 × 106
8
8, 6, 6, 6, 6, 6, 6, 6Six different numbers of MVs are specified: 3, 4, 5, 6, 7, and 8. The table shows that the requirement of the use of too few MVs will lead to additional costs when the total number of channels is the same. The cost also rises if more than 4 MVs are used, and the smaller the total number of channels is, the faster the increase of cost is. When the total number of channels is fewer than or equal to 42, the cost rises at a significantly higher rate than in the other cases. Therefore, the total number of channels should be set to greater than 42 to make it easier to achieve lower costs. However, increasing the number of channels does not significantly reduce the overall cost once the total number of channels exceeds 48. The number of MVs has a greater impact on the cost. These findings suggest that the optimization of the number of MVs is the dominant factor affecting the value of the objective function, while the number of channels is the secondary influencing factor. These observations can provide a reasonable basis for practical design.Next, a set of well distribution data from the field and a set of data randomly generated by MATLAB are used as examples to analyze how the well distribution density affects the optimization result. In the first of these examples, the coordinates of 40 wells (Table1; W69∼108) that are densely distributed in an actual production field are considered for optimization. The number of MVs m is treated as an optimization variable, and each MV is specified to have 16 channels; in the optimization result obtained, four MVs are used. It can be seen from Figure 12(a) that the optimization result is reasonable.Figure 12
Optimization results for 40 wells: (a) 40 wells in a real production field; (b) 40 wells generated randomly by MATLAB.
(a)
(b)In the second example, the positions of the 40 wells are randomly generated using the rand function (Table3; W1∼W40), and the number of MVs is constrained to a range. The upper limit on the number of channels is uniformly set to 16 for each MV, and the problem is solved using the proposed algorithm. The well connections of the six MVs are shown in Figure 12(b). These optimization results are also reasonable.


### 5.4. Stability of the Algorithm

In this section, scatter plots are presented in Figure13 to show the well and MV distribution results calculated using the proposed algorithm under the condition that the maximum number of channels per MV is constrained to 16 and the number of MVs is constrained to a certain range. The results prove that the algorithm can stably solve the MINLP problem and obtain reasonable results. Figures 13(a) and 13(b) show the result generated using the data based on W140∼W180 (Table 1), Figures 13(c) and 13(d) show the result generated using the data based on W67∼W128 (Table 1), Figure 13(e) shows the result generated using the data based on W1∼W80 (Table 1), and Figure 13(f) shows the result generated using the data based on W1∼W80 (Table 4).Figure 13
Results for various situations: (a) 41 wells, 4∼8 MVs, 16 channels per MV; (b) 41 wells, 5 MVs, 16 channels per MV; (c) 60 wells, 5 MVs, 16 channels per MV; (d) 60 wells, 6∼10 MVs, 16 channels per MV; (e) 80 wells, 5 MVs, 16 channels per MV; (f) 80 wells, 6∼10 MVs, 16 channels per MV.
(a)
(b)
(c)
(d)
(e)
(f)


### 5.5. Algorithm Performance

The scale of the computations performed in the proposed algorithm mainly depends on the size of the foresight matrix and the number of iterations. Under the assumption that all solutions are feasible, the maximum time complexity of the algorithm isomnN, where m and n denote the numbers of rows and columns, respectively, of the matrix and N is the number of iterations. If the coordinate data are given in order, the time consumption is smaller, and if they are random, the time consumption is larger, as seen by comparing the examples with 191 wells in Table 7. The random data for 191 wells are shown in Table 5. For the examples with 26 and 41 wells, for the first example of each type listed in Table 7, the identity matrix is considered the initial feasible solution, and the computational scale is then reduced through multistep iteration. By comparison with the second example in each case, it is confirmed that the multistep iteration strategy can efficiently reduce the scale of the computations performed, thus saving time. However, regardless of the type of iteration used, the same optimal solutions can be found. In the second listed examples with 35, 60, and 80 wells, there is a specified upper bound on the number of channels for each MV, and the number of MVs is constant. The computation times are increased in these cases mainly because the single-step iteration process requires more iterations to approach the optimal solution. The solution times shown in Table 7 essentially meet our expectations with regard to the algorithm complexity. Thus, it is proven that the proposed algorithm can be used to solve various problem instances within a short time.Table 7
Performance of the algorithm.
Wells
Type
Channels per MV
m (range)
f(x) (¥)
m (outcome)
Time
26
Field
All 16
2∼6
0.866 × 106
2
0.255 s
26
Field
16, 16
2
0.866 × 106
2
0.335 s
35
Field
All 16
3∼6
1.118 × 106
3
0.836 s
35
Field
8, 10, 8, 10
4
1.394 × 106
4
0.854 s
41
Field
All 16
3∼8
1.316 × 106
4
0.877 s
41
Field
All 16
4
1.316 × 106
4
1.277 s
60
Field
All 16
6∼10
1.91 × 106
6
2.554 s
60
Field
16, 8, 12, 12, 8, 10
6
1.945 × 106
6
3.418 s
80
Random
All 16
6∼10
4.330 × 106
7
5.553 s
80
Random
10, 10, 16, 16, 8, 8, 12, 12
8
4.381 × 106
8
8.336 s
191
Field
All 16
12
5.350 × 106
12
42.817 s
191
Random
All 16
12
7.382 × 106
12
78.645 sFinally, a large-scale problem with 192 wells (Table1; W1∼W192) is considered, and the result is shown in Figure 14. Here, the channel count for each valve is constrained to 16, and the number of valves is allowed to vary from 15 to 39. It can be seen that the algorithm can still successfully separate different wells into different blocks in the face of a large number of wells with multiple valves. The results show that the algorithm can still effectively solve such a large-scale problem and yield reliable results. All of the above solution processes were implemented in MATLAB on an ASUS GX501VSK laptop (Intel Core i7-7700HQ Processor, 16 GB of DDR4 RAM).Figure 14
Optimization results for a large-scale example.


## 5.1. Verification of the Validity of the Algorithm

First, a simple example is presented to demonstrate that the proposed algorithm can effectively solve the problem of interest. The associations between the wells and the MVs and the number of MVs are taken as variables, and their values are found by the algorithm via iteration; in particular, during the iterative process, multistep iteration can help the objective function value to escape from a local optimum. Let us consider a simple example for which we can make a reasonable prediction of the expected solution. Then, we can apply the proposed algorithm to solve the example and compare its solution with our expected result to prove that the algorithm can effectively solve the presented MINLP model.For all examples in this paper, we consider steel pipes with an inner diameter of 53 mm and a wall thickness of 3.5 mm as the pipes linking the wells and valves; thus, the cost coefficientwip is a constant. Moreover, the cost of an MV wjv is a known constant once the vender is decided. In this paper, the cost of an MV is treated as an average value that does not vary with the number of channels. Here, we consider an example of a field with 26 wells (Table 1; W10∼W35), for which the channel counts of the MVs and the number of MVs are not limited; therefore, the number of channels of each MV is uniformly set to 16, which is the maximum number of channels per MV that can be provided by general manufacturers.We first perform single-step iteration, during which the number of MVs will usually stay at a relatively large value, as shown in Figure7. When the identity matrix is used as the initial matrix, the number of MVs cannot be decreased below 8 via this process. However, increasing the step size can lead to further reduction in the value of the objective function, demonstrating that the reason that the objective function cannot be further reduced via single-step iteration is that the step size of each advance is too small. Thus, the objective function becomes stuck in a local optimum, and it is necessary to increase the step size to escape from that local optimum to allow the objective function to continue to decrease. The objective function decreases almost linearly during the initial 18 steps because the main factor affecting the decrease in the objective function value during the initial iterations is the number of MVs rather than the pipeline cost. However, the value of the objective function no longer decreases once it drops into a local optimum. After the 18th iteration, the step size must be increased to continue to reduce the objective function value. The single-step iteration process converges prematurely after the 18th iteration; thus, the findings demonstrate that the strategy of increasing the step size can efficiently help the objective function to escape from a local optimum.Figure 7
Comparison of two different step sizes.The result of optimization with an increase in step size is shown in Figure8. The hollow square symbols represent the positions of the MVs, the triangular symbols represent the wells connected to MV 1, and the circular symbols represent the wells connected to MV 2. For the 26 wells considered in this calculation, the algorithm completely divides them into two groups, and the solution satisfies all the constraints, which is obviously consistent with our expectations. This example proves that the proposed algorithm can be applied to solve the problem of interest and obtain reasonable optimization results.Figure 8
Results of optimization with an increase in step size.


## 5.2. Convergence of the Algorithm

Usually, we can judge whether an algorithm is valid according to whether the value of the objective function converges during the iterative process. If the algorithm is divergent, no matter how many times it is iterated, the value of the objective function will always fluctuate and will not approach a constant. In contrast, convergence of the algorithm means that the objective function value will tend toward stability after some numbers of iterations. In summary, convergence is a necessary condition for the algorithm to be valid.Next, various examples are presented to analyze the convergence of the algorithm. First, we consider a field with 35 wells (Table2; W1∼W35) as an example. MVs with four different channel counts are specified: N1up=14,N2up=12,N3up=10,andN4up=8. Thus, the number of channels of each MV is fixed in advance, the limitations on the numbers of pipelines that are allowed to be connected to different valves are different, and the number of MVs is also given, i.e., m = 4.Figure9 shows the iterative process for these 35 wells, which includes single-step iteration, double-step iteration, and redistribution of different kinds of valves. The iterative process still converges, and the optimal solution is obtained after the 21st iteration. After 14 iterations, the redistribution of different valve types is necessary to allow the objective function value to continue to decrease. Finally, a difference of 0 between two consecutive steps indicates that the final solution has been obtained.Figure 9
Iterative process for 35 wells under a fixed number of MVs with fixed channel counts.As seen from Figure9, the objective function value gradually tends toward stability and converges as iteration continues. The result is shown in Figure 10. Red circular symbols, blue triangular symbols, purple diamond symbols, and yellow triangular symbols are used to represent the individual wells connected to MVs 1, 2, 3, and 4, respectively, and square symbols are used to represent the MVs. The initial channel counts of the MVs corresponding to these different well groups are 14, 12, 10, and 8, respectively, while the final channel counts are 12, 10, 14, and 8, in sequence. In the case of tight constraints, the minimum value of the objective function can be obtained by using the proposed optimization algorithm.Figure 10
Optimization results for 35 wells.In addition to the example above, an example with 40 wells (Table1; W69∼108) and various channel counts is considered to prove the convergence of the algorithm; the results are shown in Figures 11(a) and 11(b). It can be seen from the figures that, during the initial iterations, the objective function value decreases at a relatively fast rate. Then, the rate of decrease slowly declines as iteration proceeds. The algorithm eventually converges when vN−1min−vNmin=0.Figure 11
Convergence examples under different constraints: (a) 6 valves, 40∼48 channels; (b) 8 valves, 40∼48 channels.
(a)
(b)Again, the value of the objective function gradually converges and tends toward stability. In Figure11(a), the iterative processes with a total of 46 and 48 channels are identical, indicating that once the number of channels exceeds a certain threshold, a further increase in the number of channels has no effect on the iterative process or the result. Additionally, although the iterative process for an initial feasible solution with 40 channels is different from the other 4 cases, the final result is not significantly different from the others, indicating that the initial solution does not have a great influence on the result.The iterative processes for total channel counts of 46 and 48 are also fundamentally similar. Once the total number of channels exceeds 46, the iterative processes are basically the same. The results show that the algorithm presented above can stably converge. Therefore, these examples demonstrate that even under different constraints and in different scenarios, the algorithm can drive the objective function value toward stable convergence.


## 5.3. Parameter Analysis

In general, the designer of a pipeline system should consider the possibility of reserving channels for new wells in the future along with other factors when specifying the channel counts of the MVs and the number of MVs. However, different combinations of MVs with different numbers of channels might affect the results; therefore, we will analyze the effects of different channel counts and different numbers of MVs on the results of our algorithm. We consider a field with 40 wells (Table1; W69∼108) as our example here, varying the constraints on the numbers of channels and valves to obtain different results. The results are shown in Table 6, from which we can see the influence of these different constraint parameters on the cost of the pipeline system.Table 6
Costs with various numbers of channels and MVs.
Total channels
Cost (¥)
Number of MVs
Individual channel counts
40
1.424 × 106
3
16, 16, 8
42
1.396 × 106
3
16, 16, 10
44
1.396 × 106
3
16, 14, 14
46
1.396 × 106
3
16, 16, 14
48
1.396 × 106
3
16, 16, 16
50
1.304 × 106
3
16, 16, 18
40
1.306 × 106
4
10, 10, 10, 10
42
1.287 × 106
4
12, 10, 10, 10
44
1.287 × 106
4
12, 12, 10, 10
46
1.287 × 106
4
12, 12, 12, 10
48
1.287 × 106
4
12, 12, 12, 12
50
1.287 × 106
4
14, 12, 12, 12
40
1.647 × 106
5
8, 8, 8, 8, 8
42
1.520 × 106
5
8, 8, 8, 8, 10
44
1.426 × 106
5
8, 8, 8, 10, 10
46
1.426 × 106
5
8, 8, 10, 10, 10
48
1.426 × 106
5
8, 8, 10, 10, 12
50
1.426 × 106
5
8, 8, 10, 12, 12
40
1.641 × 106
6
7, 6, 7, 7, 7, 6
42
1.619 × 106
6
8, 8, 6, 6, 6, 8
44
1.615 × 106
6
8, 8, 8, 8, 6, 6
46
1.592 × 106
6
8, 8, 8, 10, 6, 6
48
1.592 × 106
6
10, 8, 8, 10, 6, 6
50
1.592 × 106
6
10, 8, 8, 8, 8, 8
40
1.893 × 106
7
6, 6, 6, 6, 6, 4, 6
42
1.870 × 106
7
6, 6, 6, 6, 6, 6, 6
44
1.775 × 106
7
6, 8, 6, 6, 6, 6, 6
46
1.766 × 106
7
8, 8, 6, 6, 6, 6, 6
48
1.766 × 106
7
6, 8, 6, 8, 8, 6, 6
50
1.766 × 106
7
8, 8, 6, 8, 6, 6, 6
40
1.957 × 106
8
5, 5, 5, 5, 5, 5, 5, 5
42
1.950 × 106
8
6, 6, 5, 5, 5, 5, 5, 5
44
1.941 × 106
8
6, 6, 5, 5, 6, 6, 5, 5
46
1.941 × 106
8
6, 6, 5, 5, 6, 6, 6, 6
48
1.941 × 106
8
6, 6, 5, 5, 6, 6, 8, 6
50
1.919 × 106
8
8, 6, 6, 6, 6, 6, 6, 6Six different numbers of MVs are specified: 3, 4, 5, 6, 7, and 8. The table shows that the requirement of the use of too few MVs will lead to additional costs when the total number of channels is the same. The cost also rises if more than 4 MVs are used, and the smaller the total number of channels is, the faster the increase of cost is. When the total number of channels is fewer than or equal to 42, the cost rises at a significantly higher rate than in the other cases. Therefore, the total number of channels should be set to greater than 42 to make it easier to achieve lower costs. However, increasing the number of channels does not significantly reduce the overall cost once the total number of channels exceeds 48. The number of MVs has a greater impact on the cost. These findings suggest that the optimization of the number of MVs is the dominant factor affecting the value of the objective function, while the number of channels is the secondary influencing factor. These observations can provide a reasonable basis for practical design.Next, a set of well distribution data from the field and a set of data randomly generated by MATLAB are used as examples to analyze how the well distribution density affects the optimization result. In the first of these examples, the coordinates of 40 wells (Table1; W69∼108) that are densely distributed in an actual production field are considered for optimization. The number of MVs m is treated as an optimization variable, and each MV is specified to have 16 channels; in the optimization result obtained, four MVs are used. It can be seen from Figure 12(a) that the optimization result is reasonable.Figure 12
Optimization results for 40 wells: (a) 40 wells in a real production field; (b) 40 wells generated randomly by MATLAB.
(a)
(b)In the second example, the positions of the 40 wells are randomly generated using the rand function (Table3; W1∼W40), and the number of MVs is constrained to a range. The upper limit on the number of channels is uniformly set to 16 for each MV, and the problem is solved using the proposed algorithm. The well connections of the six MVs are shown in Figure 12(b). These optimization results are also reasonable.


## 5.4. Stability of the Algorithm

In this section, scatter plots are presented in Figure13 to show the well and MV distribution results calculated using the proposed algorithm under the condition that the maximum number of channels per MV is constrained to 16 and the number of MVs is constrained to a certain range. The results prove that the algorithm can stably solve the MINLP problem and obtain reasonable results. Figures 13(a) and 13(b) show the result generated using the data based on W140∼W180 (Table 1), Figures 13(c) and 13(d) show the result generated using the data based on W67∼W128 (Table 1), Figure 13(e) shows the result generated using the data based on W1∼W80 (Table 1), and Figure 13(f) shows the result generated using the data based on W1∼W80 (Table 4).Figure 13
Results for various situations: (a) 41 wells, 4∼8 MVs, 16 channels per MV; (b) 41 wells, 5 MVs, 16 channels per MV; (c) 60 wells, 5 MVs, 16 channels per MV; (d) 60 wells, 6∼10 MVs, 16 channels per MV; (e) 80 wells, 5 MVs, 16 channels per MV; (f) 80 wells, 6∼10 MVs, 16 channels per MV.
(a)
(b)
(c)
(d)
(e)
(f)


## 5.5. Algorithm Performance

The scale of the computations performed in the proposed algorithm mainly depends on the size of the foresight matrix and the number of iterations. Under the assumption that all solutions are feasible, the maximum time complexity of the algorithm isomnN, where m and n denote the numbers of rows and columns, respectively, of the matrix and N is the number of iterations. If the coordinate data are given in order, the time consumption is smaller, and if they are random, the time consumption is larger, as seen by comparing the examples with 191 wells in Table 7. The random data for 191 wells are shown in Table 5. For the examples with 26 and 41 wells, for the first example of each type listed in Table 7, the identity matrix is considered the initial feasible solution, and the computational scale is then reduced through multistep iteration. By comparison with the second example in each case, it is confirmed that the multistep iteration strategy can efficiently reduce the scale of the computations performed, thus saving time. However, regardless of the type of iteration used, the same optimal solutions can be found. In the second listed examples with 35, 60, and 80 wells, there is a specified upper bound on the number of channels for each MV, and the number of MVs is constant. The computation times are increased in these cases mainly because the single-step iteration process requires more iterations to approach the optimal solution. The solution times shown in Table 7 essentially meet our expectations with regard to the algorithm complexity. Thus, it is proven that the proposed algorithm can be used to solve various problem instances within a short time.Table 7
Performance of the algorithm.
Wells
Type
Channels per MV
m (range)
f(x) (¥)
m (outcome)
Time
26
Field
All 16
2∼6
0.866 × 106
2
0.255 s
26
Field
16, 16
2
0.866 × 106
2
0.335 s
35
Field
All 16
3∼6
1.118 × 106
3
0.836 s
35
Field
8, 10, 8, 10
4
1.394 × 106
4
0.854 s
41
Field
All 16
3∼8
1.316 × 106
4
0.877 s
41
Field
All 16
4
1.316 × 106
4
1.277 s
60
Field
All 16
6∼10
1.91 × 106
6
2.554 s
60
Field
16, 8, 12, 12, 8, 10
6
1.945 × 106
6
3.418 s
80
Random
All 16
6∼10
4.330 × 106
7
5.553 s
80
Random
10, 10, 16, 16, 8, 8, 12, 12
8
4.381 × 106
8
8.336 s
191
Field
All 16
12
5.350 × 106
12
42.817 s
191
Random
All 16
12
7.382 × 106
12
78.645 sFinally, a large-scale problem with 192 wells (Table1; W1∼W192) is considered, and the result is shown in Figure 14. Here, the channel count for each valve is constrained to 16, and the number of valves is allowed to vary from 15 to 39. It can be seen that the algorithm can still successfully separate different wells into different blocks in the face of a large number of wells with multiple valves. The results show that the algorithm can still effectively solve such a large-scale problem and yield reliable results. All of the above solution processes were implemented in MATLAB on an ASUS GX501VSK laptop (Intel Core i7-7700HQ Processor, 16 GB of DDR4 RAM).Figure 14
Optimization results for a large-scale example.


## 6. Conclusion

This paper proposes the replacement of traditional metering processes with a process based on MVs, thus simplifying and automating the necessary operations and reducing the investment in redundant valve instrumentation and other equipment, the amount of space occupied by this equipment, and the cost of the station. A generalized MINLP model is presented, and a corresponding heuristic algorithm for MV optimization design is proposed for solving the presented MINLP problem. The numbers of channels of the MVs and the number of MVs are treated as special constraints in the MINLP problem.The coordinates of actual production well sites and randomly generated data are used as examples for calculation and analysis. From the analysis of these examples, it can be seen that the behavior of the objective function value in the proposed algorithm is strictly convergent. A large collection of optimization results demonstrates that the algorithm shows stable convergence and good reliability. The analysis results show that the proposed heuristic algorithm can quickly and effectively solve the presented MINLP model for the design of a pipeline system with MVs. In this paper, we do not consider the optimization of the MV process itself in an integrated gathering system; this problem will be addressed in the future.



---

*Source: 7576093-2019-09-10.xml*