use NNLPF or NNLPG in fortran IMSL

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 · ‎06-11-2012

IMSL cannot and does not define what is and what is not legal Fortran. The Fortran standard does, and the compiler claiming standard conformance implements that definition.

The CASE-selector-block construct in Fortran allows the selector to be an initialization expression or even something such as

CASE (1, 2, 7, 10:17, 23)

You are not even required to use a SELECT construct. You can use a matrix, use binary search, or whatever you please, as long as your code does what IMSL wants it to.

Abac1 · ‎06-12-2012

Hi mecej4:Thank you for the ideas.

I tried to put my concerns into an example here. The example is adapted from example 2.8.1, pp.15-17 of Floudas's book "A Collection of Test Problems for Constrained Global Optimization Algorithms". You may find it from google book as well. (http://books.google.com.sg/books?id=zRIfcNQ0V04C&pg=PA5&lpg=PA5&dq=a+collection+of+ch2+Quadratic+Programming+test+problems&source=bl&ots=1Hh-VYhFGj&sig=cU0J89l9e_NnG11ZnDqZavyMgyo&hl=en&sa=X&ei=LxLYT8iFFsTZrQe9ueDmDw&ved=0CEQQ6AEwAA#v=onepage&q=a%20collection%20of%20ch2%20Quadratic%20Programming%20test%20problems&f=false)

I put them into the attached files called "NonOptEx1". There are a pdf file (http://software.intel.com/file/44390) and Fortran90 source file (http://software.intel.com/file/44391) inside. I have two problems with this example:

1. I cannot get the indicated solution

2. If possible, a better way coding it in Fortran should be: [min f(x)=AX^2+BX, subjected to CX=D].

In the attached source code, the problem is programmed straightforward and not good enough, especially the way of handling the constraints. I am looking for simple ways to formulate the constraints in a matrix form. Could you help to look into my attached code to make it more concise? Thanks a lot.

huang_l_ · ‎10-31-2014

i am a new learner about IMSL and i also have a similar question here.

No.1 i write my FCN subroutine as following:

SUBROUTINE FCN1 (X, IACT, RESULT, IERR)
INTEGER IACT
REAL(KIND(1E0)) X(*), RESULT
LOGICAL IERR
!
SELECT CASE (IACT)
CASE(0)
RESULT = (X(1)-2.0E0)**2 + (X(2)-1.0E0)**2
CASE(1)
RESULT = (X(1)-1)**2*(X(1))**2
!CASE(2)
!RESULT = -(X(1)**2)/4.0E0 - X(2)**2 + 1.0E0
END SELECT
RETURN
END

and i call the NNLPF in my main program like this:

   EXTERNAL FCN1
   INTEGER::M, ME
   REAL(KIND(1E0)) FVALUE, X(5), XGUESS(5), XLB(5), XUB(5)
    M=1
    ME=1
    XLB = -HUGE(X(1))
    XUB = HUGE(X(1))
    
    CALL NNLPF (FCN1, M, ME, IBTYPE, XLB, XUB, X,IPRINT=1)
    !DO i =1,size(X)
    !    write(*,'(F4.0)')X(i)
    !END DO
    
    CALL WRRRN ('The solution is', X)

as you see,my object function is "min f=(x1-2)^2+(x2-1)^2" and i only have one equality constraint:(x1-1)^2*x1^2=0; the optimal solution i think is x1=1.0,x2=1.0,and the minimun f=1.0.

but the indicated result is x1=0.0,x2=1.0,min f =4.0,as show in the attachments.

i don't know why?

No.2 if i don't know the number of equality constraints and uneqaulity constraints, and even the size of variable X, or some other linear constraints bounds ,when i write my FCN subroutine. it just decided by input files read from the main program .then how can write my FCN subroutine? in other words, how can i pass the optimization dimension and the constraints bounds parameter to my FCN?

thanks a lot!

mecej4 · ‎11-03-2014

The majority of nonlinear constrained optimization solvers are local optimizers. Their success often depends on your giving them good starting values and well-behaved functions in the objective and constraint functions. In your case, you seem to seek a global optimum, yet you gave no starting guess for X(1:2). The IMSL manual states "A good initial starting point is very problem specific and should be provided by the calling program whenever possible. See optional argument XGUESS.". You can also help move the optimizer towards the region of the more desirable solution by specifying additional constraints. For example, try x₁ >= 0.001.

What was the reason for specifying the constraint as [(x₁- 1) x₂]² = 0 instead of simply (x₁- 1) x₂ = 0 ?

huang_l_ · ‎11-04-2014

mecej4 wrote:

The majority of nonlinear constrained optimization solvers are local optimizers. Their success often depends on your giving them good starting values and well-behaved functions in the objective and constraint functions. In your case, you seem to seek a global optimum, yet you gave no starting guess for X(1:2). The IMSL manual states "A good initial starting point is very problem specific and should be provided by the calling program whenever possible. See optional argument XGUESS.". You can also help move the optimizer towards the region of the more desirable solution by specifying additional constraints. For example, try x₁ >= 0.001.

What was the reason for specifying the constraint as [(x₁- 1) x₂]² = 0 instead of simply [(x₁- 1) x₂] = 0 ?

okay, your answer is very helpful! but what about the question NO.2

"if i don't know the number of equality constraints and uneqaulity constraints, and even the size of variable X, or some other linear constraints bounds ,when i write my FCN subroutine. it just decided by input files read from the main program .then how can write my FCN subroutine? in other words, how can i pass the optimization dimension and the constraints bounds parameter to my FCN?"

thanks very much!

mecej4 · ‎11-04-2014

In situations where the argument list (specified by IMSL for the user function FCN) is not adequate, you can supplement it by using module variables (or, in the old days, using COMMON blocks). Be aware, however, that hiding information from the optimization routine can often cause failure to obtain a solution.

Here is an example, which I constructed by extending the example given in the IMSL manual for NNLPF.

    MODULE sharedvars
       IMPLICIT NONE
       INTEGER :: n, m = 2, me = 1, ibtype = 0
    END

    PROGRAM xnnlpf
       USE nnlpf_int
       USE wrrrn_int
       USE sharedvars

       IMPLICIT NONE
!
       REAL fvalue
       REAL, ALLOCATABLE :: x(:), xguess(:), xlb(:), xub(:)
       INTEGER i
       EXTERNAL fcn
!
       READ (*, *) n
       ALLOCATE (x(n), xguess(n), xlb(n), xub(n))
       xlb = -100
       xub = 100
       DO i = 1, n
          xguess(i) = 0.1*i
       END DO
!
       CALL nnlpf(fcn, m, me, ibtype, xlb, xub, x)
!
       CALL wrrrn('The solution is', x)
       DEALLOCATE (x, xguess, xlb, xub)
    END

    SUBROUTINE fcn(x, iact, result, ierr)
       USE sharedvars
       IMPLICIT NONE
       INTEGER iact, i
       REAL (kind(1E0)) x(n), & ! N is taken from SHAREDVARS        
          result
       LOGICAL ierr
!
       SELECT CASE (iact)
       CASE (0)
          result = 0
          DO i = 1, n
             result = result + (x(i)+i-3)**2
          END DO
       CASE (1)
          result = x(1) - 2.0E0*x(2) + 1.0E0
       CASE (2)
          result = -(x(1)**2)/4.0E0 - x(2)**2 + 1.0E0
       END SELECT
       RETURN
    END

Compile and run it for various values on n, such as n = 2, n=3, n=4, etc.

Andrew_Smith · ‎11-04-2014

If you contain the FCN routine inside your program it gives scope to program variables inside FCN and also provides an explicit interface to FNC so that its argument declarations can be checked against the requirements. Now we dont need the EXTERNAL statement. Note that you did not declare or set IBTYPE.

PROGRAM TEST
   use NNLPF_INT
   INTEGER IBTYPE, M, ME
   REAL(KIND(1E0)) FVALUE, X(5), XGUESS(5), XLB(5), XUB(5)
   M=1
   ME=1
   XLB = -HUGE(X(1))
   XUB = HUGE(X(1))
	     
   CALL NNLPF (FCN1, M, ME, IBTYPE, XLB, XUB, X,IPRINT=1)
   !DO i =1,size(X)
   !    write(*,'(F4.0)')X(i)
   !END DO
	     
   CALL WRRRN ('The solution is', X)
CONTAINS
   SUBROUTINE FCN1(X, IACT, RESULT, IERR)
      INTEGER, intent(in) :: IACT
      REAL(KIND(1E0)), INTENT(IN) :: X(*)
      REAL(KIND(1E0)), INTENT(OUT) :: RESULT
      LOGICAL, INTENT(OUT) :: IERR
      
      WRITE(*), XLB ! Check we can see variables in scope of our container
      SELECT CASE (IACT)
      CASE(0)
         RESULT = (X(1)-2.0E0)**2 + (X(2)-1.0E0)**2
      CASE(1)
         RESULT = (X(1)-1)**2*(X(1))**2
      !CASE(2)
      !RESULT = -(X(1)**2)/4.0E0 - X(2)**2 + 1.0E0
      END SELECT
   END SUBROUTINE
END PROGRAM

mecej4 · ‎11-05-2014

Note that you did not declare or set IBTYPE

It is declared and initialized to zero in the module.

And that is the problem with overuse of modules -- one tends to forget all the stuff put into the module, and it is ever so convenient to say "USE ...." when we no longer remember all the things that are now available. Thank goodness for "USE mymod, only: fork" and "USE mymod, except=>kitchen_sink".

huang_l_ · ‎11-06-2014

mecej4 wrote:

In situations where the argument list (specified by IMSL for the user function FCN) is not adequate, you can supplement it by using module variables (or, in the old days, using COMMON blocks). Be aware, however, that hiding information from the optimization routine can often cause failure to obtain a solution.

Here is an example, which I constructed by extending the example given in the IMSL manual for NNLPF.

    MODULE sharedvars
       IMPLICIT NONE
       INTEGER :: n, m = 2, me = 1, ibtype = 0
    END

    PROGRAM xnnlpf
       USE nnlpf_int
       USE wrrrn_int
       USE sharedvars

       IMPLICIT NONE
!
       REAL fvalue
       REAL, ALLOCATABLE :: x(:), xguess(:), xlb(:), xub(:)
       INTEGER i
       EXTERNAL fcn
!
       READ (*, *) n
       ALLOCATE (x(n), xguess(n), xlb(n), xub(n))
       xlb = -100
       xub = 100
       DO i = 1, n
          xguess(i) = 0.1*i
       END DO
!
       CALL nnlpf(fcn, m, me, ibtype, xlb, xub, x)
!
       CALL wrrrn('The solution is', x)
       DEALLOCATE (x, xguess, xlb, xub)
    END

    SUBROUTINE fcn(x, iact, result, ierr)
       USE sharedvars
       IMPLICIT NONE
       INTEGER iact, i
       REAL (kind(1E0)) x(n), & ! N is taken from SHAREDVARS        
          result
       LOGICAL ierr
!
       SELECT CASE (iact)
       CASE (0)
          result = 0
          DO i = 1, n
             result = result + (x(i)+i-3)**2
          END DO
       CASE (1)
          result = x(1) - 2.0E0*x(2) + 1.0E0
       CASE (2)
          result = -(x(1)**2)/4.0E0 - x(2)**2 + 1.0E0
       END SELECT
       RETURN
    END

Compile and run it for various values on n, such as n = 2, n=3, n=4, e

tc.

okay,I have passed the optimization parameter to my FCN as you do . but i have another question here: if i have hundreds of constraints ,for example,500 equalities constraints and 500 unequalities ,with more than 5000 variables,how can i write the FCN? i don't think the "SELECT CASE" structure is a feasible way to do this. so i try the program using the "IF -ELSEIF -END IF" statements like below:

!  shuc.f90 
!
!  FUNCTIONS:
!  shuc - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: shuc
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************
MODULE my_data_type!自定义数据类型，包括机组数据，负荷数据
  TYPE:: unit_data
      CHARACTER (len=8) :: unit_name
      CHARACTER (len=4) :: unit_zone
      REAL :: zone_trans_pmax
      REAL::unit_pmax,unit_pmin
      CHARACTER (len=1) :: unit_type
      
      TYPE(unit_data),POINTER::next_unit_data
  END TYPE  
  TYPE:: load_data
      CHARACTER (len=4) :: load_zone
      
      REAL :: peak_load,valley_load,zone_trans_pmax

      TYPE(load_data),POINTER::next_load_data
  END TYPE  
  
  
  
  TYPE:: variable


      REAL::x

      
      TYPE(variable),POINTER::next_variable
  END TYPE  
  
  
  
  
  
  
  
  
  
  CONTAINS
  SUBROUTINE unit_data_add(pos,head,unit_name,unit_zone,zone_trans_pmax,unit_pmax, unit_pmin,unit_type,unit_data_num)
      
      INTEGER::unit_data_num
      CHARACTER (len=8) :: unit_name
      CHARACTER (len=4) :: unit_zone
      REAL ::zone_trans_pmax
      REAL::unit_pmax,unit_pmin
      CHARACTER (len=1) :: unit_type
      TYPE(unit_data),POINTER::temp,pos,head
      ALLOCATE(temp)
      temp%unit_name=unit_name
      temp%unit_zone=unit_zone
      temp%zone_trans_pmax=zone_trans_pmax
      temp%unit_pmax=unit_pmax
      temp%unit_pmin=unit_pmin
      temp%unit_type=unit_type
      NULLIFY(temp%next_unit_data)

      IF (.NOT.ASSOCIATED(POS))THEN
        ALLOCATE(pos)
        pos%unit_name=unit_name
        pos%unit_zone=unit_zone
        pos%zone_trans_pmax=zone_trans_pmax
        pos%unit_pmax=unit_pmax
        pos%unit_pmin=unit_pmin
        pos%unit_type=unit_type
        
        head=>pos
        
        NULLIFY(pos%next_unit_data)
      ELSE
        pos%next_unit_data=>temp
        pos=>temp
        IF( unit_data_num==1) THEN
          head%next_unit_data=>pos
        END IF
        NULLIFY(pos%next_unit_data)
      

        

        
        

      
      
      END IF
!      DEALLOCATE( temp)
      unit_data_num=unit_data_num+1

  END  SUBROUTINE unit_data_add
  
  
  SUBROUTINE load_data_add(pos,head,load_zone,peak_load,valley_load,zone_trans_pmax,load_data_num)
      
      INTEGER::load_data_num

      CHARACTER (len=4) :: load_zone
      REAL ::zone_trans_pmax
      REAL::peak_load,valley_load

      TYPE(load_data),POINTER::temp,pos,head
      ALLOCATE(temp)
      temp%load_zone=load_zone
      temp%peak_load=peak_load
      temp%zone_trans_pmax=zone_trans_pmax
      temp%valley_load=valley_load

      NULLIFY(temp%next_load_data)

      IF (.NOT.ASSOCIATED(POS))THEN
        ALLOCATE(pos)
        pos%load_zone=load_zone
        pos%peak_load=peak_load
        pos%zone_trans_pmax=zone_trans_pmax
        pos%valley_load=valley_load

        
        head=>pos
        
        NULLIFY(pos%next_load_data)
      ELSE
        pos%next_load_data=>temp
        pos=>temp
        IF( load_data_num==1) THEN
          head%next_load_data=>pos
        END IF
        NULLIFY(pos%next_load_data)
      

        

        
        

      
      
      END IF
!      DEALLOCATE( temp)
      load_data_num=load_data_num+1

  END  SUBROUTINE load_data_add
  
    SUBROUTINE variable_add(pos,head,variable_num)
      
      INTEGER::variable_num
     


      TYPE(variable),POINTER::temp,pos,head
      ALLOCATE(temp)
      temp%x=1.

      NULLIFY(temp%next_variable)

      IF (.NOT.ASSOCIATED(POS))THEN
        ALLOCATE(pos)
        pos%x=1.

        
        head=>pos
        
        NULLIFY(pos%next_variable)
      ELSE
        pos%next_variable=>temp
        pos=>temp
        IF( variable_num==1) THEN
          head%next_variable=>pos
        END IF
        NULLIFY(pos%next_variable)
      

        

        
        

      
      
      END IF
!      DEALLOCATE( temp)
      variable_num=variable_num+1

  END  SUBROUTINE variable_add
  
  

  
  
  
  
  
  
  
END MODULE
  

MODULE  read_data
USE my_data_type
CONTAINS
SUBROUTINE read_optimization_para(peak_dc,valley_dc,peak_ac,valley_ac,head_unit,unit_data_num,head_load,load_data_num)
    CHARACTER (len=8) :: unit_name
    CHARACTER (len=4) :: unit_zone
    REAL :: zone_trans_pmax
    REAL::unit_pmax,unit_pmin
    CHARACTER (len=1) :: unit_type
    REAL::peak_dc,valley_dc,peak_ac,valley_ac
    INTEGER::i,j
    INTEGER::unit_data_num
    CHARACTER(len=1000)::data_str
    TYPE(unit_data),POINTER::head_unit,tail_unit,ptr_unit
    TYPE(unit_data)::current_unit
    LOGICAL::DATA_READING=.FALSE.
    INTEGER::result_string_scan=0
    
    INTEGER::load_data_num
    
    CHARACTER (len=4) :: load_zone
    
    REAL::peak_load,valley_load,peak_total_load,valley_total_load
    
    TYPE(load_data),POINTER::head_load,tail_load,ptr_load
    INTEGER::size_of_variable
    NULLIFY(tail_load)
    NULLIFY(head_load)
    NULLIFY(head_unit)
    NULLIFY(tail_unit)
    !NULLIFY(head_variable)
    !NULLIFY(tail_variable)
    
    OPEN(UNIT=20,FILE='上海机组小开机.txt',ACTION='READ')
    load_data_num=0
    unit_data_num=0
    DO !读入机组数据
        READ(20,'(A)') data_str
        result_string_scan=INDEX(data_str,'end_of_unit_data')
        IF (result_string_scan/=0)THEN
            DATA_READING=.FALSE.
            EXIT
        END IF
        IF(data_str(1:1)=='.')THEN
            CYCLE
        ELSE
            IF (DATA_READING==.TRUE.)THEN
                READ(data_str,*)unit_name,unit_zone,zone_trans_pmax,unit_pmax,unit_pmin,unit_type
                CALL unit_data_add(tail_unit,head_unit,unit_name,unit_zone,zone_trans_pmax,unit_pmax, unit_pmin,unit_type,unit_data_num)
                
            ELSE
                result_string_scan=INDEX(data_str,'begin_unit_data')
                IF (result_string_scan/=0)THEN
                    DATA_READING=.TRUE.
                END IF
                !result_string_scan=INDEX(data_str,'end_of_unit_data')
                !IF (result_string_scan/=0)THEN
                !    DATA_READING=.FALSE.
                !    EXIT
                !END IF
                
            END IF
            
        END IF
        
        
        
        
    
    END DO
    peak_total_load=0
    DO !读入分区负荷数据
        READ(20,'(A)') data_str
        
        result_string_scan=INDEX(data_str,'end_of_load_data')
        IF (result_string_scan/=0)THEN
            DATA_READING=.FALSE.
            EXIT
        END IF
        
        IF(data_str(1:1)=='.')THEN
            CYCLE
        ELSE
            IF (DATA_READING==.TRUE.)THEN
                READ(data_str,*)load_zone,peak_load,valley_load,zone_trans_pmax
                CALL load_data_add(tail_load,head_load,load_zone,peak_load,valley_load,zone_trans_pmax,load_data_num)
                peak_total_load=peak_total_load+peak_load
            ELSE
                result_string_scan=INDEX(data_str,'begin_load_data')
                IF (result_string_scan/=0)THEN
                    DATA_READING=.TRUE.
                END IF
                !result_string_scan=INDEX(data_str,'end_of_load_data')
                !IF (result_string_scan/=0)THEN
                !    DATA_READING=.FALSE.
                !    EXIT
                !END IF
                
            END IF
            
        END IF
        
        
        
        
    
    END DO
    
    
    DO !读入高峰低谷区外来电数据
        READ(20,'(A)') data_str
        
        result_string_scan=INDEX(data_str,'end_of_power_injection_data')
        IF (result_string_scan/=0)THEN
            DATA_READING=.FALSE.
            EXIT
        END IF
        
        IF(data_str(1:1)=='.')THEN
            CYCLE
        ELSE
            IF (DATA_READING==.TRUE.)THEN
                READ(data_str,*)peak_dc,valley_dc,peak_ac,valley_ac
                
                
            ELSE
                result_string_scan=INDEX(data_str,'begin_power_injection_data')
                IF (result_string_scan/=0)THEN
                    DATA_READING=.TRUE.
                END IF
                !result_string_scan=INDEX(data_str,'end_of_load_data')
                !IF (result_string_scan/=0)THEN
                !    DATA_READING=.FALSE.
                !    EXIT
                !END IF
                
            END IF
            
        END IF
        
        
        
        
    
    END DO
    
    CLOSE(20)
    
END SUBROUTINE read_optimization_para
END MODULE
    

    
    
    
    program shuc
    include 'link_fnl_shared.h'  
    include 'link_fnl_static.h' 
    use lin_eig_self_int
    use lin_sol_svd_int
    USE WRRRN_INT
    USE NNLPF_INT
    USE my_data_type
    USE READ_DATA

    
    
    
    
    
    
    
    IMPLICIT NONE
    INTEGER::IBTYPE=0
    INTEGER::M, ME
    CHARACTER (len=8) :: unit_name
    CHARACTER (len=4) :: unit_zone
    !REAL :: zone_trans_pmax
    !REAL::unit_pmax,unit_pmin
    !CHARACTER (len=1) :: unit_type
    INTEGER::variable_num=0
    REAL::k=1
    
    
    
    
    
    REAL::peak_dc,valley_dc,peak_ac,valley_ac!高峰低谷时段区外直流、交流来电
    
    INTEGER,PARAMETER::N=4
    !REAL::A(N,N)
    !REAL::D(N)
    INTEGER::i,j
    INTEGER::unit_data_num
    !CHARACTER(len=1000)::data_str
    
    TYPE(unit_data),POINTER::head_unit,ptr_unit
    
    !TYPE(unit_data)::current_unit
    !LOGICAL::DATA_READING=.FALSE.
    !INTEGER::result_string_scan=0
    
    INTEGER::load_data_num
    
    !CHARACTER (len=4) :: load_zone
    
    !REAL::peak_load,valley_load,peak_total_load,valley_total_load
    
    TYPE(load_data),POINTER::head_load,ptr_load
    
    
    !TYPE(variable),POINTER::head_variable,tail_variable
    
    INTEGER::size_of_variable
    REAL(KIND(1E0)) FVALUE, X(300), XGUESS(300), XLB(300), XUB(300)
    EXTERNAL FCN
    NULLIFY(ptr_load)
    NULLIFY(ptr_unit)

    
    CALL read_optimization_para(peak_dc,valley_dc,peak_ac,valley_ac,head_unit,unit_data_num,head_load,load_data_num)

    
    
    
    !ME=2!等式约束两个，分别为整数约束和高峰电力平衡约束
    !M=2*load_data_num+ME!总约束为：分区数*2（高峰、低谷分区主变不过载）+ME
    
    !!定义变量列表
    !DO i =1,3*unit_data_num
    !   CALL variable_add(tail_variable,head_variable,variable_num)
    !END DO
    
    !size_of_variable=3*unit_data_num
    !
    M=2
    ME=1
    ptr_unit=>head_unit
    DO i=1,unit_data_num
        XLB(i)=ptr_unit%unit_pmin
        XUB(i)=ptr_unit%unit_pmax
        ptr_unit=>ptr_unit%next_unit_data
        XGUESS(i)=1
        XGUESS(i+unit_data_num)=XUB(i)!迭代初始方式为高峰、低谷全开机满发
        XGUESS(i+2*unit_data_num)=XUB(i)
    END DO
    NULLIFY(ptr_unit)
    !XLB = -HUGE(X(1))
    !XUB = HUGE(X(1))
    
    CALL NNLPF (FCN, M, ME, IBTYPE, XLB, XUB, X,IPRINT=1)
    !DO i =1,size(X)
    !    write(*,'(F4.0)')X(i)
    !END DO
    
    CALL WRRRN ('The solution is', X)
    pause

    
 



    end program shuc
    
    
    
    
SUBROUTINE FCN (X, IACT, RESULT, IERR)
USE READ_DATA
USE my_data_type
INTEGER IACT
REAL(KIND(1E0)) X(*), RESULT
LOGICAL IERR



INTEGER::variable_num

REAL::peak_dc,valley_dc,peak_ac,valley_ac!高峰低谷时段区外直流、交流来电
INTEGER::i,j
INTEGER::unit_data_num=0

    
TYPE(unit_data),POINTER::head_unit,ptr_unit
    


    
INTEGER::load_data_num=0



REAL::peak_total_load,valley_total_load
    
TYPE(load_data),POINTER::head_load,ptr_load
    

TYPE(variable),POINTER::head_variable
    
NULLIFY(ptr_load)
NULLIFY(ptr_unit)
    
CALL read_optimization_para(peak_dc,valley_dc,peak_ac,valley_ac,head_unit,unit_data_num,head_load,load_data_num)
peak_total_load=0
valley_total_load=0
ptr_load=>head_load
DO i =1,load_data_num
    peak_total_load=peak_total_load+ptr_load%peak_load
    valley_total_load=valley_total_load+ptr_load%valley_load
    ptr_load=>ptr_load%next_load_data
END DO
    



variable_num=3*unit_data_num
WRITE(*,*) peak_total_load,peak_ac,peak_dc


IF(IACT==0)THEN!目标函数，低谷本地调峰缺口最小
    RESULT=0
    DO i=1,unit_data_num
        RESULT =RESULT+ X(i)*X(i+unit_data_num*2)
    END DO
    RESULT=RESULT-valley_total_load-valley_ac-valley_dc
ELSEIF(IACT==1)THEN!高峰电力平衡约束
    RESULT=0
    DO i=1,unit_data_num
        RESULT =RESULT+ X(i)*X(i+unit_data_num)
    END DO
    RESULT=RESULT-peak_total_load-peak_ac-peak_dc
    
    
ELSEIF(IACT==2)THEN!机组开机状态0-1约束
    RESULT=0
    DO i=1,unit_data_num
        RESULT=RESULT+X(i)**2*(1-X(i))**2
    END DO
ELSEIF(IACT<=3+load_data_num)THEN !高峰主变过载约束
    IF (IACT==3)THEN
        ptr_load=>head_load
        
    END IF
    RESULT=ptr_load%zone_trans_pmax-ptr_load%peak_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+unit_data_num)*X(i)
        END IF
        ptr_unit=>ptr_unit%next_unit_data
    END DO
    
    ptr_load=>ptr_load%next_load_data
    
ELSEIF(IACT<=3+2*load_data_num)THEN !低谷主变过载约束

    IF (IACT==3)THEN
        ptr_load=>head_load
        
    END IF
    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
    
    ptr_load=>ptr_load%next_load_data
END IF


RETURN
END    
    
    
    
    
    
    
    
i don't know if it will work,and there are so many "can't find or open PDB files" when linking ,and the .exe just exit without any results,shown in the following picture.

could you help me again?  thanks very much!!

Steven_L_Intel1 · ‎11-06-2014

Ignore the "can't find or open PDB file" messages when you start debugging (not linking). They are not important. I also suggest that if you want to include longer sources, please attach them as .f90 files rather than inserting into the post.

mecej4 · ‎11-06-2014

Your code compiles and links with no errors, but to run it the data file is needed.

Whether you use SELECT constructs or not is up to you. If you can group your constraints into, say, three blocks, with all the constraints in each block being of similar form, you could do something similar to

if (iact.gt.0.and.iact.le.500) then

......

else if(iact .lt. 1000) then

result = x(iact) - x(iact+400)

else

.....

endif

huang_l_ · ‎11-12-2014

mecej4 wrote:

Your code compiles and links with no errors, but to run it the data file is needed.

Whether you use SELECT constructs or not is up to you. If you can group your constraints into, say, three blocks, with all the constraints in each block being of similar form, you could do something similar to

if (iact.gt.0.and.iact.le.500) then

......

else if(iact .lt. 1000) then

result = x(iact) - x(iact+400)

else

.....

endif

i have uploaded my .f90 files and the .txt input files needed to run the program in the attachments.

but it can't get any result in the output cmd window.

i also try the "CALL NNLPF (FCN, M, ME, IBTYPE, XLB, XUB, X,XGUESS,IPRINT=1)" , where i offered my XGUESS as an initial guess of the solution, which i think is helpful to get the right answer,but i get an err as following:

error #6285: There is no matching specific subroutine for this generic subroutine call. [NNLPF]

can you help me? thanks a lot!

mecej4 · ‎11-12-2014

There are several errors in your program. You have declared several arrays that should be of size N (= the number of variables in the optimization problem) with dimension (300). The IMSL optimization routine obtains N from the size of the array passed to it, so it thinks that N = 300, whereas you want N = 21.

The input bounds must satisfy xlb < xub. You have several variables for which the bounds given satisfy xlb = xub. Do you want those elements of x() to be fixed at the bound values? If so, reformulate the problem by separating variable variables and fixed "variables".

To pass an optional argument in Fortran 90, you have to follow certain rules. To pass XGUESS, you need to write ",XGUESS=XGUESS", using the keyword name.

You should not use the == or .EQ. operators to compare logical expressions. Use .EQV. instead.

Your FCN subroutine is rather long and complicated. Have you checked whether all the components of the objective and constraints are calculated correctly? If not, I suggest that you replace the call to NNLPF with several calls "CALL FCN(xguess,Iact,fvalue,ierr)" with different values of Iact, remove all references to IMSL, and run with runtime checks turned on for uninitialized variables, array bound errors and so forth. I have a suspicion that your linked list of pointers is such that some uninitialized pointers are being dereferenced.

The routine also calls "read_optimization_para" at its beginning, which means that the data file gets read and data structures get initialized hundreds of times. Is this really necessary?

huang_l_ · ‎11-12-2014

mecej4 wrote:

There are several errors in your program. You have declared several arrays that should be of size N (= the number of variables in the optimization problem) with dimension (300). The IMSL optimization routine obtains N from the size of the array passed to it, so it thinks that N = 300, whereas you want N = 21.

The input bounds must satisfy xlb < xub. You have several variables for which the bounds given satisfy xlb = xub. Do you want those elements of x() to be fixed at the bound values? If so, reformulate the problem by separating variable variables and fixed "variables".

To pass an optional argument in Fortran 90, you have to follow certain rules. To pass XGUESS, you need to write ",XGUESS=XGUESS", using the keyword name.

You should not use the == or .EQ. operators to compare logical expressions. Use .EQV. instead.

Your FCN subroutine is rather long and complicated. Have you checked whether all the components of the objective and constraints are calculated correctly? If not, I suggest that you replace the call to NNLPF with several calls "CALL FCN(xguess,Iact,fvalue,ierr)" with different values of Iact, remove all references to IMSL, and run with runtime checks turned on for uninitialized variables, array bound errors and so forth. I have a suspicion that your linked list of pointers is such that some uninitialized pointers are being dereferenced.

The routine also calls "read_optimization_para" at its beginning, which means that the data file gets read and data structures get initialized hundreds of times. Is this really necessary?

i have modified my codes and it can run to an result now,even though the result is not what i think is correct. i guess the possible reason is that the NNLPF routine is not desighed to find global solution.

actually, in my equality contraints below:

ELSEIF(IACT==2)THEN!the 0-1constraints
    RESULT=0
    DO i=1,unit_data_num
        RESULT=RESULT+X(i)**2*(1-X(i))**2*9999
    END DO

i just want my X(1,len(X)/3)={0,1}, that is binaray numbers. but i can't find the routine that can solve optimization problem containing binaray variables in the help document. so can you tell me which routine to use to solve this problem?

thanks very much!

huang_l_ · ‎11-13-2014

supplimentary comments: my equality and unequlity constraints consists of X_i*X_{i+unit_num}items,where X_i (i>=1 and i<=unit_num )are binary variables(0-1) and X_{i+unit_num}items(_{unit_num} is a constant) are real numbers. jusu like the following constraints :

X₁*X_{1+unit_num +} X₂*X_{2+unit_num}+.....+ X_{unit_num}*X_{unit_num +unit_num}=constant₁

X₁*X_{1+unit_num +} X₂*X_{2+unit_num}>=constant₂

X₃*X_{3+unit_num +} X₄*X_{4+unit_num +}X₅*X_{5+unit_num}>=>=constant₃

X₆*X_{6+unit_num +} X₇*X_{7+unit_num +}X₈*X_{8+unit_num}+X₉*X_{9+unit_num}+X₁₀*X_{10+unit_num}>=constant₄

_.........

is there any routine to solve this problem in the IMSL package?

thanks!

mecej4 · ‎11-13-2014

If you want to do pseudo-boolean variable optimization and/or global optimization, you will have to look elsewhere for suitable algorithms. Consult Prof. Mittelmann's site at http://plato.asu.edu/guide.html .

huang_l_ · ‎11-13-2014

mecej4 wrote:

If you want to do pseudo-boolean variable optimization and/or global optimization, you will have to look elsewhere for suitable algorithms. Consult Prof. Mittelmann's site at http://plato.asu.edu/guide.html .

thanks very much!

huang_l_ · ‎11-16-2014

huang l. wrote:

Quote:

mecej4 wrote:

There are several errors in your program. You have declared several arrays that should be of size N (= the number of variables in the optimization problem) with dimension (300). The IMSL optimization routine obtains N from the size of the array passed to it, so it thinks that N = 300, whereas you want N = 21.

The input bounds must satisfy xlb < xub. You have several variables for which the bounds given satisfy xlb = xub. Do you want those elements of x() to be fixed at the bound values? If so, reformulate the problem by separating variable variables and fixed "variables".

To pass an optional argument in Fortran 90, you have to follow certain rules. To pass XGUESS, you need to write ",XGUESS=XGUESS", using the keyword name.

You should not use the == or .EQ. operators to compare logical expressions. Use .EQV. instead.

Your FCN subroutine is rather long and complicated. Have you checked whether all the components of the objective and constraints are calculated correctly? If not, I suggest that you replace the call to NNLPF with several calls "CALL FCN(xguess,Iact,fvalue,ierr)" with different values of Iact, remove all references to IMSL, and run with runtime checks turned on for uninitialized variables, array bound errors and so forth. I have a suspicion that your linked list of pointers is such that some uninitialized pointers are being dereferenced.

The routine also calls "read_optimization_para" at its beginning, which means that the data file gets read and data structures get initialized hundreds of times. Is this really necessary?

i have modified my codes and it can run to an result now,even though the result is not what i think is correct. i guess the possible reason is that the NNLPF routine is not desighed to find global solution.

actually, in my equality contraints below:
ELSEIF(IACT==2)THEN!the 0-1constraints
    RESULT=0
    DO i=1,unit_data_num
        RESULT=RESULT+X(i)**2*(1-X(i))**2*9999
    END DO
i just want my X(1,len(X)/3)={0,1}, that is binaray numbers. but i can't find the routine that can solve optimization problem containing binaray variables in the help document. so can you tell me which routine to use to solve this problem?

thanks very much!

i have another question here that why the indicated result after iteration doesn't meet the 1st equality constraints as shown below?

ELSEIF(IACT==1)THEN!高峰电力平衡约束
    RESULT=0
    DO i=1,unit_data_num
        local_p_peak =local_p_peak+ X(i)*X(i+unit_data_num)
        RESULT =RESULT+ X(i)*X(i+unit_data_num)
    END DO
    RESULT=RESULT-peak_total_load+peak_ac+peak_dc

thank you again!

mecej4 · ‎11-16-2014

why the indicated result after iteration doesn't meet the 1st equality constraints as shown below

The IMSL routine gave up on the calculation and returned without any iterations. In fact, it did not even evaluate the object function once before giving up. See this part of the output:

TERMINATION REASON:
STEPSIZESELECTION: NO ACCEPTABLE STEPSIZE IN [SIGSM,SIGLA]
EVALUATIONS OF F 0
EVALUATIONS OF GRAD F 0
EVALUATIONS OF CONSTRAINTS 16535
EVALUATIONS OF GRADS OF CONSTRAINTS 0

I would not be brave enough to try such a large problem without having good starting data. You have set the bounds as (-1,2) for many variables, and the starting value as 1. Are these reasonable? Can you do better? I am not going to be of much help in this part because I have no knowledge of the application area. Do the variables have units (metres, pascals, volts, etc.)?

mecej4 · ‎11-18-2014

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".