use NNLPF or NNLPG in fortran IMSL - Page 2

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.

huang_l_ · ‎11-19-2014

mecej4 wrote:

Do you have a version of the problem with fewer variables and constraints? If so, you may be able to solve that problem first, and use the results from it to generate good starting values for the larger problem.

After thinking and looking at your problem over a few hours, I think that I have found a way to solve it. Here are my thoughts.

You have N = 3 n variables, where n is what you call unit_data_num in the program.Let the unknowns be called X = {u, v, w} where u = X(1:n), v = X(n+1:2n), w = X(2n+1:N). Then, your objective function is

F = u^Tw + a

and the general constraints are of the form

some(u.*v) >= b, some other(u.*w) >= c

In order to force u to be an array of binary variables you have added the artificial constraint

norm {u.*(1 - u)} = 0

The binary variables are not only going to cause standard continuous optimization methods to have a lot of trouble, but they are also unnecessary!

Use new variables {y, z}, where y = u.*v and z = u.*w, where ".*" stands for element-by-element multiplication, as in Matlab. With these variables, you no longer need the additional constraint, and you have a third fewer variables. The problem that you have implemented in the program may now be recast as a linear programming (LP) problem involving continuous variables only, and could be solved efficiently using an LP solver.

I have attached a file containing the solution that I found (using a modified version of your program, which I have attached also). In the solution file, if a line in the file ends with an 'l', that variable is at the lower bound; if a 'u', it is at the upper bound. Otherwise, the variable is free (bounds not active). You may interpret those variables that are zero to be entities such as "power plant P is on but its output is 0 or power plant P is off but if it had been on then the output would have been 200 MW but not delivered to consumer/customer Q".

Thank you very much for help!

The optimal solution after iteration in your modified program is the "maximun scenario" in the real electrical power grid world, which means that all the power plants is on during a typical day. So in your result,each element in X(1:88) is greater than 0. But in my optimization model, I just want to find an solution that contain zeros in X(1:88) ,as you newly defined. It will get an better objective value in FCN, and the solution really works in the running power system currently. But it has an implicit constraint in your X(2*unit_data_num):if X(i)==0,the X(i+num) has to be equal to 0. which means if output of a coal fueled plant in day time is 0, then it has to be off line, and further more the X(i+unit_data_num) ,the output of a coal fueled plant at night, has to be 0 due to its off status.That is why I have to use binary variables.

huang_l_ · ‎11-19-2014

sorry,there are some errors in my comments above, it should be :

"But it has an implicit constraint in your X(2*unit_data_num):if X(i)==0,the X(i+unit_data_num) has to be equal to 0"

mecej4 · ‎11-19-2014

Thanks for the explanations; now I understand the problem better, and need to spend some time to think about the implications.

I have an immediate question. You stated that if a plant is off during the day, it must also be off at night. Is the converse also true? I can imagine that if a plant is coal-fired and off during the day it cannot be on at night but, if you have a number of peak load plants such as gas-turbine or diesel plants that can be quickly brought on line, it should be acceptable that many or all X(i+unit_data_num) = 0 even though the corresponding X(i) > 0.

If so, the following solution may be acceptable. An interior optimizer such as NNLPF will find a local optimum close to the true optimum, whereas a linear programming algorithm may do better at finding the exact "vertex" that gives an optimum. For your problem, using the LP solver SLPRS in IMSL, the solution is found to be

Plant   1 day on  night off    2535.00
Plant   8 day on  night off    3773.30
Plant  28 day on  night off     929.10
Plant  37 day on  night off     193.50
Plant  55 day on  night off    1082.85
Plant  68 day on  night off    1717.25
Plant  77 day on  night off    1200.70
Plant  81 day on  night off    3947.70

Please comment about this "solution", in which X[1, 8, 28, 37, 55, 68, 77, 81] are the only non-zero values in the optimum solution. I realize that this may not be a meaningful engineering solution because the objective function value, which is Σ z_i- 7072, is negative, but this is a consequence of the current mathematical formulation.

There seems to be another set of equations in your original formulation that I cannot understand. You have two instances of code such as this:

    RESULT=ptr_load%zone_trans_pmax-ptr_load%valley_load
    ptr_unit=>head_unit
    DO i =1,unit_data_num
       IF (ptr_unit%unit_zone==ptr_load%load_zone)THEN
           RESULT=RESULT+X(i+2*unit_data_num)*X(i)
       END IF
       ptr_unit=>ptr_unit%next_unit_data
    END DO

Because of the IF..ENDIF block inside the loop, for IACT=12 to 15 and 25 to 28, RESULT does not receive any contributions from any of the X variables. This is probably because the last four load zone names in the data file do not match any of the zones in the units section. Mathematically, this situation creates a constraint of the type 0 >= - RESULT. You may need to add code to omit such constraints, which cannot be satisfied. Since RESULT is not zero even if the DO loop is not executed, how is the load met with no plant running?

huang_l_ · ‎11-20-2014

mecej4 wrote:

Thanks for the explanations; now I understand the problem better, and need to spend some time to think about the implications.

I have an immediate question. You stated that if a plant is off during the day, it must also be off at night. Is the converse also true? I can imagine that if a plant is coal-fired and off during the day it cannot be on at night but, if you have a number of peak load plants such as gas-turbine or diesel plants that can be quickly brought on line, it should be acceptable that many or all X(i+unit_data_num) = 0 even though the corresponding X(i) > 0.

If so, the following solution may be acceptable. An interior optimizer such as NNLPF will find a local optimum close to the true optimum, whereas a linear programming algorithm may do better at finding the exact "vertex" that gives an optimum. For your problem, using the LP solver SLPRS in IMSL, the solution is found to be
Plant   1 day on  night off    2535.00
Plant   8 day on  night off    3773.30
Plant  28 day on  night off     929.10
Plant  37 day on  night off     193.50
Plant  55 day on  night off    1082.85
Plant  68 day on  night off    1717.25
Plant  77 day on  night off    1200.70
Plant  81 day on  night off    3947.70
Please comment about this "solution", in which X[1, 8, 28, 37, 55, 68, 77, 81] are the only non-zero values in the optimum solution. I realize that this may not be a meaningful engineering solution because the objective function value, which is Σ z_i- 7072, is negative, but this is a consequence of the current mathematical formulation.

There seems to be another set of equations in your original formulation that I cannot understand. You have two instances of code such as this:
    RESULT=ptr_load%zone_trans_pmax-ptr_load%valley_load
    ptr_unit=>head_unit
    DO i =1,unit_data_num
       IF (ptr_unit%unit_zone==ptr_load%load_zone)THEN
           RESULT=RESULT+X(i+2*unit_data_num)*X(i)
       END IF
       ptr_unit=>ptr_unit%next_unit_data
    END DO
Because of the IF..ENDIF block inside the loop, for IACT=12 to 15 and 25 to 28, RESULT does not receive any contributions from any of the X variables. How do you explain this? Since RESULT is not zero even if the DO loop is not executed, how is the load met with no plant running?

yes, you are right! in the real power system ,if a plant is off during the day, it must also be off at night and vice versa, that is why i use extral binary variables to represent the status of power plants in the optimization model.

i am sorry to say that the optimal solution abtained from SLPRS subroutine is not a meaningful solution. actually, the unit data type in the "shuc.txt" has an "unit_type" attribute , which is defined in the MODULE my_data_type. and "C" means coal mine plant;"Z" means enterprise owned plants and "G" means gas turbine plant. the "C" type plant meets the constraint of "if a plant is off during the day, it must also be off at night and vice versa"; the "Z" type has to be online for all the time;and the "G" type plant has no on/off constraints, in other words ,the output of "G" type plant is just an continous variable between its lower and upper bound. but to make things simple, i don't take these 3 types of plants into consideration, and i just make an assumption that we only have the "C" type plants.

in the IF..ENDIF block you mentioned above, i just want to consider the following constraints:

ptr_load%zone_trans_pmax-ptr_load%peak_load+sum(X(i+unit_data_num)*X(i))>=0 (for any plants meeting:ptr_unit%unit_zone==ptr_load%load_zone )

and :

ptr_load%zone_trans_pmax-ptr_load%valley_load+sum(X(i+2*unit_data_num)*X(i))(for any plants meeting:ptr_unit%unit_zone==ptr_load%load_zone )

i can't understand that why "for IACT=12 to 15 and 25 to 28, RESULT does not receive any contributions from any of the X variables"?

can you help me more? thanks!!

mecej4 · ‎11-20-2014

i can't understand that why "for IACT=12 to 15 and 25 to 28, RESULT does not receive any contributions from any of the X variables"

This happens because the data is such that the test

IF (ptr_unit%unit_zone==ptr_load%load_zone)

is not satisfied for those values of IACT.

Please examine your data file (by the way, the link that you provided in #13 to the attachment shuc.txt does not work anymore). The copy that I have shows the zone name for the last four load data as I have listed below (sorry, I cannot read Chinese and the text editors that I have used may have corrupted the text). The first column shows the auto-detect rendering of the editor, the second column shows the GB18030 version, and the third column is from an online transliterator.

ÈýÁÖ 三林 sān lín
¾²°² 静安 Jìng ān 
ÐÂÓà 新余 Xīn yú
Á·ÌÁ 练塘 liàn táng

None of these names match any of the zone names in the unit data section.

Similarly, the last eight lines in the unit data section show the zone name "Ö÷Íø" or "主网" or "zhǔ wǎng", which does not appear anywhere in the load data section, but this may be fine.

huang_l_ · ‎11-21-2014

mecej4 wrote:

Quote:

i can't understand that why "for IACT=12 to 15 and 25 to 28, RESULT does not receive any contributions from any of the X variables"

This happens because the data is such that the test
IF (ptr_unit%unit_zone==ptr_load%load_zone)
is not satisfied for those values of IACT.

Please examine your data file (by the way, the link that you provided in #13 to the attachment shuc.txt does not work anymore). The copy that I have shows the zone name for the last four load data as I have listed below (sorry, I cannot read Chinese and the text editors that I have used may have corrupted the text). The first column shows the auto-detect rendering of the editor, the second column shows the GB18030 version, and the third column is from an online transliterator.
ÈýÁÖ 三林 sān lín
¾²°² 静安 Jìng ān 
ÐÂÓà 新余 Xīn yú
Á·ÌÁ 练塘 liàn táng
None of these names match any of the zone names in the unit data section.

Similarly, the last eight lines in the unit data section show the zone name "Ö÷Íø" or "主网" or "zhǔ wǎng", which does not appear anywhere in the load data section, but this may be fine.

i have submitted a "shuc.txt" file in english. can you further explain why "is not satisfied for those values of IACT."? thanks!

mecej4 · ‎11-21-2014

When I click on the link to shuc_eng.txt in #27, I get a "Page not found" notice. I still have your original shuc.txt, so let us use that for the following explanation. I thought that I had described why you have some improper constraints, but let us try again.

Consider that FCN is called with IACT=12. The code that gets executed is

        ptr_load=>head_load
        DO i=2,IACT-2
            ptr_load=>ptr_load%next_load_data
        END DO

The DO loop is executed with i = 2,3,...,10, i.e., nine times, and ptr_load points to the tenth load, which is given on line-103 of the data file. On that line, we have

三林 851.55 459.837 1150

The load-zone name, 三林 (san-lin), does not match any of the zone names in the units section of the data file, in lines 3 to 90. Therefore, the statements within the IF block inside the DO loop are not executed even once for this value of IACT.

    DO i =1,unit_data_num
        IF (ptr_unit%unit_zone==ptr_load%load_zone)THEN
            local_p_peak=local_p_peak+X(i+unit_data_num)*X(i)
            RESULT=RESULT+X(i+unit_data_num)*X(i)
        END IF
        ptr_unit=>ptr_unit%next_unit_data
    END DO

Therefore, RESULT simply has the initial value assigned to it before the DO loop is started, and is independent of the X values.

huang_l_ · ‎11-21-2014

mecej4 wrote:

When I click on the link to shuc_eng.txt in #27, I get a "Page not found" notice. I still have your original shuc.txt, so let us use that for the following explanation. I thought that I had described why you have some improper constraints, but let us try again.

Consider that FCN is called with IACT=12. The code that gets executed is
        ptr_load=>head_load
        DO i=2,IACT-2
            ptr_load=>ptr_load%next_load_data
        END DO
The DO loop is executed with i = 2,3,...,10, i.e., nine times, and ptr_load points to the tenth load, which is given on line-103 of the data file. On that line, we have

三林 851.55 459.837 1150

The load-zone name, 三林 (san-lin), does not match any of the zone names in the units section of the data file, in lines 3 to 90. Therefore, the statements within the IF block inside the DO loop are not executed even once for this value of IACT.
    DO i =1,unit_data_num
        IF (ptr_unit%unit_zone==ptr_load%load_zone)THEN
            local_p_peak=local_p_peak+X(i+unit_data_num)*X(i)
            RESULT=RESULT+X(i+unit_data_num)*X(i)
        END IF
        ptr_unit=>ptr_unit%next_unit_data
    END DO
Therefore, RESULT simply has the initial value assigned to it before the DO loop is started, and is independent of the X values.

you are absolutely right! there are no corresponding plant units in the last 4 load zones,which means "Therefore, RESULT simply has the initial value assigned to it before the DO loop is started, and is independent of the X values."

i have condition checking codes in my program as shown below:

    ptr_load=>head_load
    DO i=1,load_data_num
        zone_pmax=0
        ptr_unit=>head_unit
        DO j=1,unit_data_num
            IF (ptr_unit%unit_zone==ptr_load%load_zone)THEN
                zone_pmax=zone_pmax+ptr_unit%unit_pmax
            END IF
            ptr_unit=>ptr_unit%next_unit_data
        END DO
        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
            
            
        END IF
        ptr_load=>ptr_load%next_load_data
        
    END DO

there are 4 load zones without any plants in it as you find. but it dosenot matter anything if it pass the condition checking process above.

huang_l_ · ‎11-21-2014

mecej4 wrote:

When I click on the link to shuc_eng.txt in #27, I get a "Page not found" notice. I still have your original shuc.txt, so let us use that for the following explanation. I thought that I had described why you have some improper constraints, but let us try again.

Consider that FCN is called with IACT=12. The code that gets executed is
        ptr_load=>head_load
        DO i=2,IACT-2
            ptr_load=>ptr_load%next_load_data
        END DO
The DO loop is executed with i = 2,3,...,10, i.e., nine times, and ptr_load points to the tenth load, which is given on line-103 of the data file. On that line, we have

三林 851.55 459.837 1150

The load-zone name, 三林 (san-lin), does not match any of the zone names in the units section of the data file, in lines 3 to 90. Therefore, the statements within the IF block inside the DO loop are not executed even once for this value of IACT.
    DO i =1,unit_data_num
        IF (ptr_unit%unit_zone==ptr_load%load_zone)THEN
            local_p_peak=local_p_peak+X(i+unit_data_num)*X(i)
            RESULT=RESULT+X(i+unit_data_num)*X(i)
        END IF
        ptr_unit=>ptr_unit%next_unit_data
    END DO
Therefore, RESULT simply has the initial value assigned to it before the DO loop is started, and is independent of the X values.

i post the context in the "shuc_eng.txt" in the comments below and you can copy it :

begin_unit_data
.unit          zone       zone_trans_pmax         pmax          pmin       type
hushxi_3      xuha       1350                    660           330        C
hushxi_4      xuha       1350                    660           330        C
hushra_1      xuha       1350                    400           0          G
hushra_2      xuha       1350                    400           0          G
hushra_3      xuha       1350                    400           0          G
hujidi_1      xuha       1350                    11            10         Z
hujidi_2      xuha       1350                    11            10         Z
huxiba_0      yaha       1000                    150           75         Z
huxiba_3      yaha       1000                    350           175        Z
huxiba_4      yaha       1000                    350           175        Z
hubaga_1      yaha       1000                    350           175        Z
hubaga_2      yaha       1000                    350           175        Z
hushyi_1      yaha       1000                    325           162.5      C
hushyi_2      yaha       1000                    325           162.5      C
hushyi_3      yaha       1000                    325           162.5      C
hushyi_4      yaha       1000                    325           162.5      C
huzhra_1      yaha       1000                    100           0          G
huzhra_2      yaha       1000                    100           0          G
huzhra_3      yaha       1000                    100           0          G
huzhra_4      yaha       1000                    100           0          G
huzhra_5      yaha       1000                    111           0          G
huzhra_6      yaha       1000                    111           0          G
huzhbe15      yaha       1000                    115           55         G
huzhbe16      yaha       1000                    115           55         G
huwega__      yaha       1000                    50            49         Z
hugadi_1      yaha       1000                    12            11         Z
hugadi_2      yaha       1000                    15            14         Z
huwayi_1      gulu       3000                    320           160        C
huwayi_2      gulu       3000                    320           160        C
huwayi_3      gulu       3000                    320           160        C
huwayi_4      gulu       3000                    320           160        C
hugaqi_4      gulu       3000                    50            49         Z
hugaqi_5      gulu       3000                    25            24         Z
hugaqi_6      gulu       3000                    25            24         Z
hugaqi_7      gulu       3000                    25            24         Z
hugaqi_8      gulu       3000                    50            49         Z
hulira_1      yudo       1350                    400           0          G
hulira_2      yudo       1350                    400           0          G
hulira_3      yudo       1350                    400           0          G
hulira_4      yudo       1350                    400           0          G
hufera31      yudo       1350                    120           0          G
hufera41      yudo       1350                    120           0          G
hufera32      yudo       1350                    60            0          G
hufera42      yudo       1350                    60            0          G
huxidi_1      yudo       1350                    11            10         Z
huxidi_2      yudo       1350                    11            10         Z
hufera11      yaga       2600                    120           0          G
hufera12      yaga       2600                    60            0          G
hufera21      yaga       2600                    120           0          G
hufera22      yaga       2600                    60            0          G
huzhlo_1      yaga       2600                    11            10         Z
huzhlo_2      yaga       2600                    11            10         Z
hupuch_1      yaga       2600                    8             7          Z
hupuch_2      yaga       2600                    8             7          Z
hujime_0      tiwe       1350                    25            24         Z
hujime_1      tiwe       1350                    50            49         Z
hujime_2      tiwe       1350                    50            49         Z
hujime_3      tiwe       1350                    50            49         Z
hujime_4      tiwe       1350                    50            49         Z
hujime_5      tiwe       1350                    100           99        Z
hujime_6      tiwe       1350                    100           99        Z
hujime_7      tiwe       1350                    100           99        Z
hujime_8      tiwe       1350                    100           99        Z
hucare11      tiwe       1350                    242           240        Z
hucare12      tiwe       1350                    90            89         Z
hucare21      tiwe       1350                    242           241        Z
hucare22      tiwe       1350                    90            89         Z
humira_1      naqi       1900                    400           0          G
humira_2      naqi       1900                    400           0          G
huwuji_8      naqi       1900                    300           299        Z
huwuji_9      naqi       1900                    300           299        Z
huwuji11      naqi       1900                    300           150        C
huwuji12      naqi       1900                    300           150        C
hujihu_1      naqi       1900                    12            11         Z
hujihu_2      naqi       1900                    12            11         Z
hujihu_3      naqi       1900                    12            11         Z
huwuer_1      siji       2600                    620           310        C
huwuer_2      siji       2600                    620           310        C
huqidi_1      hudu       2800                    11            10         Z
huqidi_2      hudu       2800                    11            10         Z
hushca_1      zhwa       9999                    1050          525        C
hushca_2      zhwa       9999                    1050          525        C
huwasa_7      zhwa       9999                    1050          525        C
huwasa_8      zhwa       9999                    1050          525        C
huwaer_5      zhwa       9999                    950           475        C
huwaer_6      zhwa       9999                    950           475        C
husher_1      zhwa       9999                    660           330        C
husher_2      zhwa       9999                    660           330        C
end_of_unit_data
begin_load_data
.zone peak_load    valley_load    zone_trans_pmax
xuha 3885        2097.9          1350
yaha 4773.3      2863.98         1000
gulu 3929.1      2121.714        3000
yudo 1543.5      833.49          1350
yaga 1907.85      1030.239        2600
tiwe 2432.85      1313.739        1350
naqi 3617.25      1953.315        1900
siji 3800.7      2052.378        2600
hudu 2730        1474.2          2800
sali 851.55      459.837         1150
jian 806.4        435.456         2700
xiyu 612.15      330.561         1350
lita 1065.75      575.505         1350
end_of_load_data
begin_power_injection_data
.peak_dc valley_dc peak_ac valley_ac
9800    9800    6776      670
end_of_power_injection_data

mecej4 · ‎11-22-2014

Thanks for the English version of the data file, and here is some progress.

Since (i) NNLPF was not finding an initial feasible point, (ii) your code did not actually implement the logic rules for different types of units, (iii) I felt that the way you used "x(1-x)" constraints to treat binary variables could cause the optimizer problems, and (iv) your code used linked lists to store and retrieve the input data, which is painful to watch in a debugger, I did two things.

I created an AMPL model for the problem to try different models to generate an initial guess for NNLPF. Among the solvers available for AMPL, few allow 0-1 variables, and many require the Hessian to be positive definite at the initial point, or do not allow nonlinear equality constraints. However, with the Gurobi solver for AMPL, by relaxing your first constraint to a pair of inequality constraints, 0 < C1 < 8.6 , I was able to obtain a solution. This solution gave F = 2959.6 as the optimum.
Using the solution from AMPL as a starting value, I called NNLPF using a substantially modified version of your program, and obtained an improved solution, F = 2436.6. The modifications include (i) using arrays instead of linked lists, (ii) imposing the different rules for generation units properly (as I understand them), and (iii) using a different barrier function instead of the sum of squares of X(1-X) of the binary variables.

Here are the function values at the optimum:

Iact    F(X,iact) Old Iact
----   ---------- --------
 0     2436.6015   0       Objective Function
 1        0.0000   1       Equality Constraint
 2        0.0000   2       Dummy Constraint (for binary variables)
 3        7.0000   3       Physical Constraint 1
 4        5.7000   4
 5      218.1625   5
 6       28.1506   6
 7      730.1500   7
 8      206.1500   8
 9      150.1037   9
10       39.3000  10
11       92.0000  11
12        0.0001  16
13        1.5879  17
14     1688.2860  18
15      536.9575  19
16     1603.7610  20
17     1311.2610  21
18      877.6850  22
19     1167.6220  23
20     1345.8000  24

        12 to 15 and 25 to 28:  These are constraints with no X terms.

I am surprised but pleased that NNLPF was able to run to completion on a problem which it was not designed for.

Please examine the results critically and let me know how you wish to proceed.

huang_l_ · ‎11-22-2014

mecej4 wrote:

Thanks for the English version of the data file, and here is some progress.

Since (i) NNLPF was not finding an initial feasible point, (ii) your code did not actually implement the logic rules for different types of units, (iii) I felt that the way you used "x(1-x)" constraints to treat binary variables could cause the optimizer problems, and (iv) your code used linked lists to store and retrieve the input data, which is painful to watch in a debugger, I did two things.

I created an AMPL model for the problem to try different models to generate an initial guess for NNLPF. Among the solvers available for AMPL, few allow 0-1 variables, and many require the Hessian to be positive definite at the initial point, or do not allow nonlinear equality constraints. However, with the Gurobi solver for AMPL, by relaxing your first constraint to a pair of inequality constraints, 0 < C1 < 8.6 , I was able to obtain a solution. This solution gave F = 2959.6 as the optimum.

Using the solution from AMPL as a starting value, I called NNLPF using a substantially modified version of your program, and obtained an improved solution, F = 2436.6. The modifications include (i) using arrays instead of linked lists, (ii) imposing the different rules for generation units properly (as I understand them), and (iii) using a different barrier function instead of the sum of squares of X(1-X) of the binary variables.

Here are the function values at the optimum:
Iact    F(X,iact) Old Iact
----   ---------- --------
 0     2436.6015   0       Objective Function
 1        0.0000   1       Equality Constraint
 2        0.0000   2       Dummy Constraint (for binary variables)
 3        7.0000   3       Physical Constraint 1
 4        5.7000   4
 5      218.1625   5
 6       28.1506   6
 7      730.1500   7
 8      206.1500   8
 9      150.1037   9
10       39.3000  10
11       92.0000  11
12        0.0001  16
13        1.5879  17
14     1688.2860  18
15      536.9575  19
16     1603.7610  20
17     1311.2610  21
18      877.6850  22
19     1167.6220  23
20     1345.8000  24

        12 to 15 and 25 to 28:  These are constraints with no X terms.
I am surprised but pleased that NNLPF was able to run to completion on a problem which it was not designed for.

Please examine the results critically and let me know how you wish to proceed.

can you let me know the elements of your optimal solution and its engineering meaning ?so i can verify whether it is the right answer for the problem. thanks a lot!

mecej4 · ‎11-23-2014

The attached zip file contains all the files needed to produce a solution. Here are the steps to reproduce my results.

Unzip the zip file into a separate directory. Move all the *.txt files (with the exception of shuce.txt) elsewhere for comparison later.
Obtain the command line version of AMPL (http://ampl.com/dl/demo/ampl-demo-mswin.zip) if you do not already have it.
Run AMPL on the model file shuc.mod. The Gurobi solver will produce the solution in a second or two. Copy/paste the solution into "asol.txt". Remove the first five lines and the last line from 'asol.txt'. You should have 88 lines left. This solution does not satisfy all the constraints that you wish to apply, and we are going to use it only as an improved "guess" in the next step.
Compile and link sharnlpf.f90 with IMSL. Run the program. It will run the optimization three times, with successively higher multiplier values applied to the barrier term added to the objective function. It will produce the optimized solution within a minute, and write the solution into "nsol.txt".
Compile and link "evalsol.f90". Now run the resulting EXE, which will read the solution file output in Step 4 and check that the solution satisfies the constraints. The output should resemble the results displayed in #31.

The only point that will remain uncertain is whether there are other local minima than the one found by the NNLPF solver. It may be the case that the one found is acceptable.

huang_l_ · ‎11-24-2014

mecej4 wrote:

The attached zip file contains all the files needed to produce a solution. Here are the steps to reproduce my results.

Unzip the zip file into a separate directory. Move all the *.txt files (with the exception of shuce.txt) elsewhere for comparison later.

Obtain the command line version of AMPL (http://ampl.com/dl/demo/ampl-demo-mswin.zip) if you do not already have it.

Run AMPL on the model file shuc.mod. The Gurobi solver will produce the solution in a second or two. Copy/paste the solution into "asol.txt". Remove the first five lines and the last line from 'asol.txt'. You should have 88 lines left. This solution does not satisfy all the constraints that you wish to apply, and we are going to use it only as an improved "guess" in the next step.

Compile and link sharnlpf.f90 with IMSL. Run the program. It will run the optimization three times, with successively higher multiplier values applied to the barrier term added to the objective function. It will produce the optimized solution within a minute, and write the solution into "nsol.txt".

Compile and link "evalsol.f90". Now run the resulting EXE, which will read the solution file output in Step 4 and check that the solution satisfies the constraints. The output should resemble the results displayed in #31.

The only point that will remain uncertain is whether there are other local minima than the one found by the NNLPF solver. It may be the case that the one found is acceptable.

i am so sorry that i do not know AMPL at all. i can't understand "Run AMPL on the model file shuc.mod" .it seem i am required to enter something to run the "ampl.exe". but i don't know what to enter?

i have an alternaltive: can you just offer the result abtained by the AMPL program to me in .txt format,then i can deal use it as a guess solution in your program by myself.

mecej4 · ‎11-24-2014

The command to run AMPL is "ampl shuc.mod".

As I wrote in #33, the results from the AMPL run are already included: file "asol.txt", so if you don't have AMPL just keep asol.txt in the working directory, then compile and run the optimization program as I described in Step-4.

huang_l_ · ‎11-26-2014

mecej4 wrote:

The command to run AMPL is "ampl shuc.mod".

As I wrote in #33, the results from the AMPL run are already included: file "asol.txt", so if you don't have AMPL just keep asol.txt in the working directory, then compile and run the optimization program as I described in Step-4.

to make things easier first, i modified your program by deleting the unit_type consideration , and makes assumption that all units are coal fired generator.

but i think there are still some glitches in the program :

No.1:i find that the optimal result depends on the guess solution. if i change the guess solution by set some status of unit to 1 (the original value is 0)and reduce other online units output to meet the equal constraints , so i get a new guess solution, and i will get an different optimal solution after running the program. this means that the iteration solution is not the global solution we want. how can we get it?

i submitted 2 guess solutions , in "asol1.txt" and "asol2.txt" ; and the 2 corresponding results are shown in "nsol1.txt" and "nsol2.txt".

No2:in the real power system word, it is very hard to find a feasible solution manully as a guess solution meeting all the constrains, (or else we don't need to use this program). So we just want to know if there is some way we can find the optimal solution without offering a "good guess", or can we get a "good guess" by any other subroutines just as the AMPL do?

mecej4 · ‎11-26-2014

Sorry, the forum software is preventing me from seeing the .TXT files that you attached. I can, however, see the .f90 file. I hope Intel will fix this problem soon. In the meantime, let us try something different about file names. (i) zip all the files in #36 and attach the zip file, with the .ZIP ending, to a new post; (ii) rename the .TXT files as .TXT.f90 and post them as attachments.

In general, the success of solving an NLP (non linear programming) problem is heavily dependent on the quality of the initial guesses. See, for example, this section from a book: http://books.google.com/books?id=BbHVAQAAQBAJ&pg=PA10&lpg=PA10&dq=hill+climbing+multiple+peaks&source=bl&ots=HZc7CBN7Vc&sig=7b9s91W8B6yGgjhHH9jdGfjegJQ&hl=en&sa=X&ei=Ltx1VPjDPIukNqClg4gK&ved=0CFMQ6AEwBw#v=onepage&q=hill%20climbing%20multiple%20peaks&f=false .

to make things easier first, i modified your program by deleting the unit_type consideration , and makes assumption that all units are coal fired generator.

Let me warn you about such simplifications, as you may not be aware of the consequences, such as the following. In the actual application, as you have stated, gas turbine plants can be turned off at night. Consider units 25 and 26, "huzhbe15" and "huzhbe16". For these units, you have lower bounds given as 55. If you now treat these units as if they are coal units, those coal units cannot be off, even at night, without violating this lower bound. It can happen that this unreasonable and unnecessarily high lower bound can result in the problem becoming infeasible. To overcome this problem, you could use X(i)*Z(i) >= ~~Xlb(i)~~ X(i)*Xlb(i) as a constraint instead of X(i) >= Xlb(i), but then you would have a general nonlinear constraint replacing a simple lower bound constraint, and that causes all sorts of difficulties. Here is one simple test that you can do to highlight this issue. Consider a simple model in which all the units are C type, and so can be on or off, but you still apply the lower bounds on Z even to units that are off. The resulting problem is a linear programming problem in new variables U = XY and V=XZ. With the inappropriate bounds applied to U and V, this problem is infeasible, as you can find by running a LP solver on the modified problem.

Steven_L_Intel1 · ‎11-26-2014

I can see the .txt files even if not logged in to the forum, so it isn't a permissions problem.

mecej4 · ‎11-26-2014

Thanks, Steve, and I understand that you cannot fix a problem that you cannot see.

From my end, the problem is that whether I am logged in or not, in Chrome or IE, I get the "page not found" page when I click on one of the links to the .txt files. In fact, using one of those links with Cygwin/GNU wget I find

s:\IMSL\coal>wget https://software.intel.com/sites/default/files/managed/6c/47/nsol2.txt
--2014-11-26 12:24:05--  https://software.intel.com/sites/default/files/managed/6c/47/nsol2.txt
Resolving software.intel.com (software.intel.com)... 23.208.112.37
Connecting to software.intel.com (software.intel.com)|23.208.112.37|:443... connected.
HTTP request sent, awaiting response... 404 Not Found
2014-11-26 12:24:06 ERROR 404: Not Found.

Just as I can download the .f90 file from a browser, wget is able to download that file.

Steven_L_Intel1 · ‎11-26-2014

Oh - the problem is that the files don't download, not that you can't see the links. THAT I can see. I will report it to the forum maintainers when I get back to the office.

mecej4 · ‎11-26-2014

Similar problems have been reported in the MKL Forum. Ying H (Intel's) post, https://software.intel.com/en-us/forums/topic/531014#comment-1798254 , refers to a zip file at https://software.intel.com/en-us/articles/using-intel-mkl-in-your-c-program . Following the link and accepting a license form takes one to https://software.intel.com/en-us/protected-download/267270/157506/step2 , but clicking on this link brings up "page not found". Ying H (Intel) responded to the complaint by attaching the files in questions as attachments in a reply.

Steve, I am a bit confused about:

1. "I can see the .txt files even if not logged in to the forum, so it isn't a permissions problem."

and

2. "the problem is that the files don't download, not that you can't see the links"

I can see the links, but not the "linkee"s of type .txt, in the forum posts, for example, the link nsol1.txt in #36 under "Attachments list". When I place the mouse on the link, the browser shows the URL in the status box as https://software.intel.com/sites/default/files/managed/6c/47/nsol1.txt, which appears to be a reasonable link to a file. My browser settings are such that common file type such as .txt and .jpg are handled by the browser and .pdf by a browser plug-in. There is no handler for .f90, so the browser will attempt to download such a file when I click on a link to it. It seems to me that the .txt links point to files that do not exist and, if so, the forum software may have a bug in forming the URL that it includes as a link in a post after an attachment has been uploaded to the Intel server.

huang_l_ · ‎11-27-2014

mecej4 wrote:

Sorry, the forum software is preventing me from seeing the .TXT files that you attached. I can, however, see the .f90 file. I hope Intel will fix this problem soon. In the meantime, let us try something different about file names. (i) zip all the files in #36 and attach the zip file, with the .ZIP ending, to a new post; (ii) rename the .TXT files as .TXT.f90 and post them as attachments.

In general, the success of solving an NLP (non linear programming) problem is heavily dependent on the quality of the initial guesses. See, for example, this section from a book: http://books.google.com/books?id=BbHVAQAAQBAJ&pg=PA10&lpg=PA10&dq=hill+c... .

Quote:

to make things easier first, i modified your program by deleting the unit_type consideration , and makes assumption that all units are coal fired generator.

Let me warn you about such simplifications, as you may not be aware of the consequences, such as the following. In the actual application, as you have stated, gas turbine plants can be turned off at night. Consider units 25 and 26, "huzhbe15" and "huzhbe16". For these units, you have lower bounds given as 55. If you now treat these units as if they are coal units, those coal units cannot be off, even at night, without violating this lower bound. It can happen that this unreasonable and unnecessarily high lower bound can result in the problem becoming infeasible. To overcome this problem, you could use X(i)*Z(i) >= Xlb(i) as a constraint instead of X(i) >= Xlb(i), but then you would have a general nonlinear constraint replacing a simple lower bound constraint, and that causes all sorts of difficulties. Here is one simple test that you can do to highlight this issue. Consider a simple model in which all the units are C type, and so can be on or off, but you still apply the lower bounds on Z even to units that are off. The resulting problem is a linear programming problem in new variables U = XY and V=XZ. With the inappropriate bounds applied to U and V, this problem is infeasible, as you can find by running a LP solver on the modified problem.

i am afraid i can't agree with you about "Consider units 25 and 26, "huzhbe15" and "huzhbe16". For these units, you have lower bounds given as 55. If you now treat these units as if they are coal units, those coal units cannot be off, even at night, without violating this lower bound. It can happen that this unreasonable and unnecessarily high lower bound can result in the problem becoming infeasible".

actually,both huzhbe15 and huzhbe16 can be shut down during a whole day just by setting the value of X(i), i<=unit_data_num ,to 0, even if X(i+unit_data_num)>=55; the program should be able to find the units that have a 0 status.the binary variables is used to represent the on/off status of a unit , no matter its lower/upper bounds. and if i use X(i)*Z(i) >= Xlb(i) as a constraint instead of X(i) >= Xlb(i), than it really implys that the unit is always online , which is not true.

i submitted a .rar file containing the 4 files yuo can't see before. please help me again.