use NNLPF or NNLPG in fortran IMSL - Page 3

Abac1 · ‎06-11-2012

In NNLPF, the constrained nonlinear optimization problem is cast into the manner of "SELECT CASE ". CASE (0) corresponds to the objective function evaluation, while CASE (1:ME) for equality constraints and CASE(ME+1:M) for unequality contraints.

If ME/M is a small number, say less than 5, we can easily write the constraints one by one in the same way as the documentation examples (Ch 8 Optimization of the IMSL documentation). However, What if ME/M is considerably large as 1000, are there any better ways to write these constraints? For example, is it possible to write the constraints via matrix? Thanks in advance.

mecej4 · ‎11-27-2014

and if i use X(I)*Z(I) >= Xlb(I) as a constraint instead of X(I) >= Xlb(I), than it really implies that the unit is always online , which is not true.

I made a slip in #37, which I have corrected, and meant to state X(I)*Z(I) >= X(I)*Xlb(I) instead. Let us analyze this formulation. X(I) is, for coal and gas stations, either 0 or 1. For such stations, if X(I)=0, the condition is satisfied regardless of the values of Z(I) and Xlb(I). For the same stations, if X(i)=1, the condition is equivalent to Z(I) >= Xlb(I). Thus, the correct bound constraints are applied (in the engineering sense; in the mathematical model, the product of X*Z is not a variable, so we cannot call the constraint a bound on a mathematical variable).

I have modified the AMPL model to implement these modified constraints, and a solution is obtained within a few seconds. Surprisingly, the results also satisfy the requirement that gas units can be operated only during the day or only during the night or both times, whichever serves to satisfy the load constraints and minimize the objective. There was no need to run NNLPF on this modified problem, since AMPL does not report any initial infeasibility. The results are given below. Please try to answer two questions that puzzle me but which, with your experience, you may explain easily:

There are several overloaded transmission zones. In your original program, you had a loop to check such overloads, but had commented out the two WRITE statements in the loop. Unfortunately, I could not translate the text of those output strings from Chinese.
The objective value is negative. To me, this implies an error in the problem formulation. You have not directly stated the engineering description of the objective function, but to me it appears to be wasted/unusable power output at night, which is to be minimized. A negative value would mean that, far from having surplus production, there is a deficit. If so, the problem formulation needs to be repaired.

I have downloaded and looked at the contents of the .rar file of #37. Both solutions (nsol1.txt and nsol2.txt) are valid ~~and give the same final value~~ but give different values of the objective function. I have also observed elsewhere that different solutions may give the same values for some of the objective and constraint functions, and differ for others. It is not surprising that when you have objective functions and constraint functions that are unaffected by interchanging one or more pairs of variables. That is, if x1 and x2 are two variables in the problem, and the objective function and all of those constraints that contain x1 or x2 contain x1 + x2, any combination of x1 and x2, subject to the bounds on those variables, will yield the same set of values for the objective function and constraint functions. Which of these combinations is a calculation going to converge to? That depends on the initial guess and the optimization algorithm.

S:\IMSL\shucn>evalsoln shucx.sol
 ncoal, ngas, nent =           22          25          41

 Overloaded Transmission Zones
           2  (yaha)
           4  (yudo)
           6  (tiwe)
           7  (naqi)
           8  (siji)

 Station Outputs
  No.   Name    C_day       P_day   C_night  P_Night Type
   1  hushxi_3   1.0       660.000           330.000  C
   2  hushxi_4   1.0       660.000           330.000  C
   3  hushra_1   1.0       400.000                    G
   4  hushra_2   1.0       400.000                    G
   5  hushra_3   1.0       400.000   1.0      67.900  G
   6  hujidi_1              11.000            10.000  Z
   7  hujidi_2              11.000            10.000  Z
   8  huxiba_0             150.000            75.000  Z
   9  huxiba_3             350.000           175.000  Z
  10  huxiba_4             350.000           175.000  Z
  11  hubaga_1             350.000           175.000  Z
  12  hubaga_2             350.000           175.000  Z
  13  hushyi_1   1.0       325.000           162.500  C
  14  hushyi_2   1.0       325.000           162.500  C
  15  hushyi_3   1.0       325.000           162.500  C
  16  hushyi_4   1.0       325.000           162.500  C
  17  huzhra_1   1.0       100.000                    G
  18  huzhra_2   1.0       100.000   1.0      34.980  G
  19  huzhra_3   1.0       100.000                    G
  20  huzhra_4   1.0       100.000   1.0     100.000  G
  21  huzhra_5   1.0       111.000                    G
  22  huzhra_6   1.0       111.000                    G
  23  huzhbe15   1.0       115.000   1.0     115.000  G
  24  huzhbe16   1.0       115.000   1.0     115.000  G
  25  huwega__              50.000            49.000  Z
  26  hugadi_1              12.000            11.000  Z
  27  hugadi_2              15.000            14.000  Z
  28  huwayi_1   1.0       320.000           160.000  C
  29  huwayi_2   1.0       320.000           160.000  C
  30  huwayi_3   1.0       320.000           160.000  C
  31  huwayi_4   1.0       306.400           160.000  C
  32  hugaqi_4              50.000            49.000  Z
  33  hugaqi_5              25.000            24.000  Z
  34  hugaqi_6              25.000            24.000  Z
  35  hugaqi_7              25.000            24.000  Z
  36  hugaqi_8              50.000            49.000  Z
  37  hulira_1   1.0       400.000                    G
  38  hulira_2   1.0       400.000                    G
  39  hulira_3   1.0       400.000                    G
  40  hulira_4   1.0       400.000                    G
  41  hufera31   1.0       120.000                    G
  42  hufera41   1.0       120.000                    G
  43  hufera32   1.0        60.000                    G
  44  hufera42   1.0        60.000                    G
  45  huxidi_1              11.000            10.000  Z
  46  huxidi_2              11.000            10.000  Z
  47  hufera11   1.0       120.000                    G
  48  hufera12   1.0        60.000                    G
  49  hufera21   1.0       120.000                    G
  50  hufera22   1.0        60.000                    G
  51  huzhlo_1              11.000            10.000  Z
  52  huzhlo_2              11.000            10.000  Z
  53  hupuch_1               8.000             7.000  Z
  54  hupuch_2               8.000             7.000  Z
  55  hujime_0              25.000            24.000  Z
  56  hujime_1              50.000            49.000  Z
  57  hujime_2              50.000            49.000  Z
  58  hujime_3              50.000            49.000  Z
  59  hujime_4              50.000            49.000  Z
  60  hujime_5             100.000            99.000  Z
  61  hujime_6             100.000            99.000  Z
  62  hujime_7             100.000            99.000  Z
  63  hujime_8             100.000            99.000  Z
  64  hucare11             242.000           240.000  Z
  65  hucare12              90.000            89.000  Z
  66  hucare21             242.000           241.000  Z
  67  hucare22              90.000            89.000  Z
  68  humira_1   1.0       400.000                    G
  69  humira_2   1.0       400.000                    G
  70  huwuji_8             300.000           299.000  Z
  71  huwuji_9             300.000           299.000  Z
  72  huwuji11   1.0       300.000           150.000  C
  73  huwuji12   1.0         0.000             0.000  C
  74  hujihu_1              12.000            11.000  Z
  75  hujihu_2              12.000            11.000  Z
  76  hujihu_3              12.000            11.000  Z
  77  huwuer_1   1.0       620.000           310.000  C
  78  huwuer_2   1.0       620.000           310.000  C
  79  huqidi_1              11.000            10.000  Z
  80  huqidi_2              11.000            10.000  Z
  81  hushca_1   1.0         0.000             0.000  C
  82  hushca_2   1.0         0.000             0.000  C
  83  huwasa_7   1.0         0.000             0.000  C
  84  huwasa_8   1.0         0.000             0.000  C
  85  huwaer_5   1.0       950.000           475.000  C
  86  huwaer_6   1.0         0.000             0.000  C
  87  husher_1   1.0         0.000             0.000  C
  88  husher_2   1.0         0.000             0.000  C

Objective:

   0     -425.4340

Constraints:

   1        0.0000
   2        0.0000
   3        7.0000
   4        5.7000
   5      512.3000
   6     1788.5000
   7     1090.1500
   8      206.1500
   9       18.7500
  10       39.3000
  11       92.0000
  12       -0.0000
  13        0.0000
  14     1688.2860
  15      536.5100
  16     1603.7610
  17     1311.2610
  18      727.6850
  19     1167.6220
  20     1345.8000

huang_l_ · ‎11-28-2014

mecej4 wrote:

Quote:

and if i use X(I)*Z(I) >= Xlb(I) as a constraint instead of X(I) >= Xlb(I), than it really implies that the unit is always online , which is not true.

I made a slip in #37, which I have corrected, and meant to state X(I)*Z(I) >= X(I)*Xlb(I) instead. Let us analyze this formulation. X(I) is, for coal and gas stations, either 0 or 1. For such stations, if X(I)=0, the condition is satisfied regardless of the values of Z(I) and Xlb(I). For the same stations, if X(i)=1, the condition is equivalent to Z(I) >= Xlb(I). Thus, the correct bound constraints are applied (in the engineering sense; in the mathematical model, the product of X*Z is not a variable, so we cannot call the constraint a bound on a mathematical variable).

I have modified the AMPL model to implement these modified constraints, and a solution is obtained within a few seconds. Surprisingly, the results also satisfy the requirement that gas units can be operated only during the day or only during the night or both times, whichever serves to satisfy the load constraints and minimize the objective. There was no need to run NNLPF on this modified problem, since AMPL does not report any initial infeasibility. The results are given below. Please try to answer two questions that puzzle me but which, with your experience, you may explain easily:

There are several overloaded transmission zones. In your original program, you had a loop to check such overloads, but had commented out the two WRITE statements in the loop. Unfortunately, I could not translate the text of those output strings from Chinese.

The objective value is negative. To me, this implies an error in the problem formulation. You have not directly stated the engineering description of the objective function, but to me it appears to be wasted/unusable power output at night, which is to be minimized. A negative value would mean that, far from having surplus production, there is a deficit. If so, the problem formulation needs to be repaired.

I have downloaded and looked at the contents of the .rar file of #37. Both solutions (nsol1.txt and nsol2.txt) are valid and give the same final value of the objective function. It is not surprising that when you have objective functions and constraint functions that are unaffected by interchanging one or more pairs of variables. That is, if x1 and x2 are two variables in the problem, and the objective function and all of those constraints that contain x1 or x2 contain x1 + x2, any combination of x1 and x2, subject to the bounds on those variables, will yield the same set of values for the objective function and constraint functions. Which of these combinations is a calculation going to converge to? That depends on the initial guess and the optimization algorithm.
S:\IMSL\shucn>evalsoln shucx.sol
 ncoal, ngas, nent =           22          25          41

 Overloaded Transmission Zones
           2  (yaha)
           4  (yudo)
           6  (tiwe)
           7  (naqi)
           8  (siji)

 Station Outputs
  No.   Name    C_day       P_day   C_night  P_Night Type
   1  hushxi_3   1.0       660.000           330.000  C
   2  hushxi_4   1.0       660.000           330.000  C
   3  hushra_1   1.0       400.000                    G
   4  hushra_2   1.0       400.000                    G
   5  hushra_3   1.0       400.000   1.0      67.900  G
   6  hujidi_1              11.000            10.000  Z
   7  hujidi_2              11.000            10.000  Z
   8  huxiba_0             150.000            75.000  Z
   9  huxiba_3             350.000           175.000  Z
  10  huxiba_4             350.000           175.000  Z
  11  hubaga_1             350.000           175.000  Z
  12  hubaga_2             350.000           175.000  Z
  13  hushyi_1   1.0       325.000           162.500  C
  14  hushyi_2   1.0       325.000           162.500  C
  15  hushyi_3   1.0       325.000           162.500  C
  16  hushyi_4   1.0       325.000           162.500  C
  17  huzhra_1   1.0       100.000                    G
  18  huzhra_2   1.0       100.000   1.0      34.980  G
  19  huzhra_3   1.0       100.000                    G
  20  huzhra_4   1.0       100.000   1.0     100.000  G
  21  huzhra_5   1.0       111.000                    G
  22  huzhra_6   1.0       111.000                    G
  23  huzhbe15   1.0       115.000   1.0     115.000  G
  24  huzhbe16   1.0       115.000   1.0     115.000  G
  25  huwega__              50.000            49.000  Z
  26  hugadi_1              12.000            11.000  Z
  27  hugadi_2              15.000            14.000  Z
  28  huwayi_1   1.0       320.000           160.000  C
  29  huwayi_2   1.0       320.000           160.000  C
  30  huwayi_3   1.0       320.000           160.000  C
  31  huwayi_4   1.0       306.400           160.000  C
  32  hugaqi_4              50.000            49.000  Z
  33  hugaqi_5              25.000            24.000  Z
  34  hugaqi_6              25.000            24.000  Z
  35  hugaqi_7              25.000            24.000  Z
  36  hugaqi_8              50.000            49.000  Z
  37  hulira_1   1.0       400.000                    G
  38  hulira_2   1.0       400.000                    G
  39  hulira_3   1.0       400.000                    G
  40  hulira_4   1.0       400.000                    G
  41  hufera31   1.0       120.000                    G
  42  hufera41   1.0       120.000                    G
  43  hufera32   1.0        60.000                    G
  44  hufera42   1.0        60.000                    G
  45  huxidi_1              11.000            10.000  Z
  46  huxidi_2              11.000            10.000  Z
  47  hufera11   1.0       120.000                    G
  48  hufera12   1.0        60.000                    G
  49  hufera21   1.0       120.000                    G
  50  hufera22   1.0        60.000                    G
  51  huzhlo_1              11.000            10.000  Z
  52  huzhlo_2              11.000            10.000  Z
  53  hupuch_1               8.000             7.000  Z
  54  hupuch_2               8.000             7.000  Z
  55  hujime_0              25.000            24.000  Z
  56  hujime_1              50.000            49.000  Z
  57  hujime_2              50.000            49.000  Z
  58  hujime_3              50.000            49.000  Z
  59  hujime_4              50.000            49.000  Z
  60  hujime_5             100.000            99.000  Z
  61  hujime_6             100.000            99.000  Z
  62  hujime_7             100.000            99.000  Z
  63  hujime_8             100.000            99.000  Z
  64  hucare11             242.000           240.000  Z
  65  hucare12              90.000            89.000  Z
  66  hucare21             242.000           241.000  Z
  67  hucare22              90.000            89.000  Z
  68  humira_1   1.0       400.000                    G
  69  humira_2   1.0       400.000                    G
  70  huwuji_8             300.000           299.000  Z
  71  huwuji_9             300.000           299.000  Z
  72  huwuji11   1.0       300.000           150.000  C
  73  huwuji12   1.0         0.000             0.000  C
  74  hujihu_1              12.000            11.000  Z
  75  hujihu_2              12.000            11.000  Z
  76  hujihu_3              12.000            11.000  Z
  77  huwuer_1   1.0       620.000           310.000  C
  78  huwuer_2   1.0       620.000           310.000  C
  79  huqidi_1              11.000            10.000  Z
  80  huqidi_2              11.000            10.000  Z
  81  hushca_1   1.0         0.000             0.000  C
  82  hushca_2   1.0         0.000             0.000  C
  83  huwasa_7   1.0         0.000             0.000  C
  84  huwasa_8   1.0         0.000             0.000  C
  85  huwaer_5   1.0       950.000           475.000  C
  86  huwaer_6   1.0         0.000             0.000  C
  87  husher_1   1.0         0.000             0.000  C
  88  husher_2   1.0         0.000             0.000  C

Objective:

   0     -425.4340

Constraints:

   1        0.0000
   2        0.0000
   3        7.0000
   4        5.7000
   5      512.3000
   6     1788.5000
   7     1090.1500
   8      206.1500
   9       18.7500
  10       39.3000
  11       92.0000
  12       -0.0000
  13        0.0000
  14     1688.2860
  15      536.5100
  16     1603.7610
  17     1311.2610
  18      727.6850
  19     1167.6220
  20     1345.8000

i don't know way you say "There are several overloaded transmission zones". your solution is perfect! it is exactly what i am looking for! the constraints value from No.3 to No.20 is greater than 0 really means that no zone is overloaded. it seems we have got the right answer now, but not by IMSL subroutines....it is done by AMPL , which i don't know at all..

but i don't agree with you saying "Both solutions (nsol1.txt and nsol2.txt) are valid and give the same final value of the objective function". i see the final value is different , 2543.56 in nsol1.txt and 2312.56 in nso21.txt. the ∑X(i)*X(i+2*unit_data_num) i <=unit_data_num is 9615 and 9384 respectively.

mecej4 · ‎11-28-2014

i don't know way you say "There are several overloaded transmission zones".

In your original file shuc.f90, you had the following lines, but with the WRITE statements commented out.

        IF(ptr_load%zone_trans_pmax<ABS(ptr_load%peak_load-zone_pmax))THEN
            WRITE(*,'(A,"    分区原始数据有误,高峰分区电力平衡无解：")')ptr_load%load_zone
            WRITE(*,'("高峰分区负荷：",F6.0,"  高峰分区最大出力：",F6.0,"  分区主变限额：",F6.0)')ptr_load%peak_load,zone_pmax,ptr_load%zone_trans_pmax

I placed these lines in the newer program as a check on the input data that you gave in file shuce.txt. I realize that one should perhaps use the results from the optimizer instead of the peak capacities of the units. These WRITE statements issue a warning that if the plants do operate at peak capacity the transmission zone capacities listed may be exceeded.

but i don't agree with you saying "Both solutions (nsol1.txt and nsol2.txt) are valid and give the same final value of the objective function". i see the final value is different

The numbers that you compare are the values of expressions calculated from the optimal values of the variables after the optimization has been completed. These expressions were not used in determining the optimum. It is a mathematical fact that if, after solving the optimization problem

f(x) = minimum, such that g_e(x) = 0, g_i(x) >= 0, x_l <= x <= x_u,

you then compute another function h(x), two valid solutions of the optimization problem, i.e., solutions which give the same f_min while satisfying all the constraints, may still give different values of h(x).

huang_l_ · ‎11-28-2014

mecej4 wrote:

Quote:

i don't know way you say "There are several overloaded transmission zones".

In your original file shuc.f90, you had the following lines, but with the WRITE statements commented out.
        IF(ptr_load%zone_trans_pmax<ABS(ptr_load%peak_load-zone_pmax))THEN
            WRITE(*,'(A,"    分区原始数据有误,高峰分区电力平衡无解：")')ptr_load%load_zone
            WRITE(*,'("高峰分区负荷：",F6.0,"  高峰分区最大出力：",F6.0,"  分区主变限额：",F6.0)')ptr_load%peak_load,zone_pmax,ptr_load%zone_trans_pmax
I placed these lines in the newer program as a check on the input data that you gave in file shuce.txt. I realize that one should perhaps use the results from the optimizer instead of the peak capacities of the units. These WRITE statements issue a warning that if the plants do operate at peak capacity the transmission zone capacities listed may be exceeded.

the following statements:

WRITE(*,'(A," 分区原始数据有误,高峰分区电力平衡无解：")')ptr_load%load_zone

WRITE(*,'("高峰分区负荷：",F6.0," 高峰分区最大出力：",F6.0," 分区主变限额：",F6.0)')ptr_load%peak_load,zone_pmax,ptr_load%zone_trans_pmax

can be translated into:

WRITE(*,'(A," load data err,there are no feasible solutions to meet the peak load scenario：")')ptr_load%load_zone

WRITE(*,'("peak load data ：",F6.0," maximun unit power in the zone ：",F6.0," maximun thansformer capacity of the zone：",F6.0)')ptr_load%peak_load,zone_pmax,ptr_load%zone_trans_pmax

these statements are initial condition checking codes to aviod "bad model" that has no feasible solution .

mecej4 wrote:

Quote:

f(x) = minimum, such that g_e(x) = 0, g_i(x) >= 0, x_l <= x <= x_u,

you then compute another function h(x), two valid solutions of the optimization problem, i.e., solutions which give the same f_min while satisfying all the constraints, may still give different values of h(x).

what is h(x)?

mecej4 · ‎11-29-2014

I have made a correction to #43. The two solutions are valid, but differ. However, those results are for the simplified "coal rules for all" model.

h(x) is any function other than those used in the optimization problem. I looked briefly at the part of your code that comes after the optimizer finishes. Since these lines of code did not call the FCN subroutine I concluded that they formed h(x). Now that you have translated the WRITE statements for me, I went back to the code and I see that the lines in question duplicate the calculation done in FCN. Therefore, although my remarks about a post-optimization calculation (of any h(x)) are correct, they do not apply to the results that you posted in nsol1.txt and nsol2.txt. The language barrier is partly to blame for this "noise".

Let us now focus on the future. The modified mathematical model that I described in #43 gives results that you have stated as "perfect". So far, we do not have an implementation of that model that runs with NNLPF. In fact, I thought about creating one, but decided against doing so for the following reasons.

AMPL works nicely (as does another optimization solver package, GAMS)with the MINLP (Mixed Integer NonLinear Program) formulation, and is very convenient to use. Changes to the mathematical model are possible with slight effort.
The modified model involves bounds on products of the optimization variables and not on the variables themselves. Implementing this model in NNLPF would significantly increase the problem size, consume a significant time for rewriting code, and slow down the optimization run. Since we know that NNLPF is really for NLP problems rather than MINLP problems, the modification is probably not worthwhile.

I urge you to try either AMPL or GAMS. If you are doing the optimization for an employer, they should realize that the expense of obtaining AMPL or GAMS is easily justified by considering how much time you are going to save by using such tailored solvers instead of trying to continue with IMSL.

huang_l_ · ‎11-29-2014

mecej4 wrote:

I have made a correction to #43. The two solutions are valid, but differ. However, those results are for the simplified "coal rules for all" model.

h(x) is any function other than those used in the optimization problem. I looked briefly at the part of your code that comes after the optimizer finishes. Since these lines of code did not call the FCN subroutine I concluded that they formed h(x). Now that you have translated the WRITE statements for me, I went back to the code and I see that the lines in question duplicate the calculation done in FCN. Therefore, although my remarks about a post-optimization calculation (of any h(x)) are correct, they do not apply to the results that you posted in nsol1.txt and nsol2.txt. The language barrier is partly to blame for this "noise".

Let us now focus on the future. The modified mathematical model that I described in #43 gives results that you have stated as "perfect". So far, we do not have an implementation of that model that runs with NNLPF. In fact, I thought about creating one, but decided against doing so for the following reasons.

AMPL works nicely (as does another optimization solver package, GAMS) with the MINLP (Mixed Integer NonLinear Program) formulation, and is very convenient to use. Changes to the mathematical model are possible with slight effort.

The modified model involves bounds on products of the optimization variables and not on the variables themselves. Implementing this model in NNLPF would significantly increase the problem size, consume a significant time for rewriting code, and slow down the optimization run. Since we know that NNLPF is really for NLP problems rather than MINLP problems, the modification is probably not worthwhile.

I urge you to try either AMPL or GAMS. If you are doing the optimization for an employer, they should realize that the expense of obtaining AMPL or GAMS is easily justified by considering how much time you are going to save by using such tailored solvers instead of trying to continue with IMSL.

thank you very much! you really helped a lot!

mecej4 · ‎11-30-2014

I appreciate your comments.

Should you wish to try GAMS, the necessary files to run the optimization are inside the attached Zip file. GAMS comes with a larger selection of MINLP solvers than AMPL.

You can run GAMS on the power system model with the command

     gams.exe SHUC.gms minlp=couenne

You can try other MINLP solvers by replacing "couenne" with "sbb", etc., in this command line. You can also use the GAMS IDE, if you wish.

After GAMS has found the optimum, the results are placed into the text file "shuc.sol". You can then run the GAMS version of evalsol.f90 on that file to print out a report on the results.