From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001 From: stainer_t Date: Mon, 8 Sep 2025 13:48:49 +0200 Subject: Initial commit from Polytechnique Montreal --- Donjon/src/PLMAP1.f | 341 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 341 insertions(+) create mode 100644 Donjon/src/PLMAP1.f (limited to 'Donjon/src/PLMAP1.f') diff --git a/Donjon/src/PLMAP1.f b/Donjon/src/PLMAP1.f new file mode 100644 index 0000000..3f4d677 --- /dev/null +++ b/Donjon/src/PLMAP1.f @@ -0,0 +1,341 @@ +*DECK PLMAP1 + SUBROUTINE PLMAP1(N0,M0,APLUS,PDG,BPLUS,INPLUS,XDROIT,COUT,OBJ, + > XOBJ,IMTHD,IMPR,IERR,BINF,BSUP,SCALE,PX,RX,DELTA,BGAR) +* +*----------------------------------------------------------------------- +* +*Purpose: +* Solves a linear optimization problem with quadratic constraint using +* the method of approximation programming (MAP). +* PLMAP1 = Linear Programmation MAP1 +* +*Reference: +* R.E. Griffith and R.A. Stewart, 'A non-linear programming technique +* for the optimization of continuous processing systems', Management +* Science, Vol. 7, NO. 4, 379 (1961). +* +*Copyright: +* Copyright (C) 2002 Ecole Polytechnique de Montreal +* This library is free software; you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation; either +* version 2.1 of the License, or (at your option) any later version +* +*Author(s): +* A. Hebert and R. Chambon +* +*Parameters: input +* N0 number of control variables. +* M0 number of constraints. +* APLUS coefficient matrix for the linear constraints. +* PDG weights assigned to control variables in the quadratic +* constraint. +* BPLUS right hand sides corresponding to the coefficient matrix. +* INPLUS constraint relations (=-1 for .GE.; =0 for .EQ.; =1 for .LE.). +* XDROIT quadratic constraint radius squared. +* COUT costs of control variables. +* OBJ objective function. +* XOBJ control variables. +* IMTHD type of solution (=1: SIMPLEX/LEMKE; =3: MAP). +* IMPR print flag. +* +*Parameters: ouput +* IERR return code (=0: normal completion). +* +*Parameters: scratch +* BINF +* BSUP +* SCALE +* PX +* RX +* DELTA +* BGAR +* +*----------------------------------------------------------------------- +* + IMPLICIT DOUBLE PRECISION (A-H,O-Z) +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER N0,M0,INPLUS(M0+1),IMTHD,IMPR,IERR + DOUBLE PRECISION PDG(N0),BPLUS(M0+2),XDROIT,XOBJ(N0),BINF(N0), + > BSUP(N0),SCALE(N0),PX(N0),RX(N0),DELTA(N0),BGAR(M0+1), + > APLUS(M0+2,N0),COUT(N0),OBJ +*---- +* LOCAL VARIABLES +*---- + DOUBLE PRECISION DELF,X,DUMY,ZMAX,XVAL,SCAL,DELX,TEMP,ERR, + > CONT,EPSIR,EPS,EPSS + INTEGER ITER,ITMAX,I,J,M + CHARACTER CLNAME*6 + DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: ZMAXC,BGAR0 +*---- +* DATA STATEMENTS +*---- + DATA CLNAME /'PLMAP'/ +* + EPS = 0.0001D0 + EPSIR =EPS + ITMAX = 100 +*---- +* CONTROL-VARIABLE SCALING +*---- + DO 10 J=1,N0 + SCAL = SQRT(XDROIT/PDG(J)) + SCALE(J) = SCAL + COUT(J) = COUT(J)*SCAL +* + PDG(J) = 1.0D0 + BINF(J) = 0.0D0 + BSUP(J) = 0.0D0 +* + DO 20 I=1,M0 + APLUS(I,J) = APLUS(I,J)*SCAL + 20 CONTINUE +* + 10 CONTINUE +*---- +* PRINT TABLES AFTER SCALING OF CONTROL VARIABLES +*---- + IF(IMPR.GE.5) THEN + CALL PLNTAB(COUT,APLUS,INPLUS,BPLUS,PDG,BINF,BSUP, + > N0,M0,CLNAME//' AFTER SCALING OF CONTROL VARIABLES') + ENDIF +* + XDROIT=1.0D0 +*---- +* CONSTRAINT SCALING +*---- + ALLOCATE(ZMAXC(M0)) + DO 30 I=1,M0 + ZMAX = ABS(BPLUS(I)) +* + DO 40 J=1,N0 + ZMAX = MAX(ZMAX,ABS(APLUS(I,J))) + 40 CONTINUE + BGAR(I) = BPLUS(I)/ZMAX +* + DO 42 J=1,N0 + APLUS(I,J) = APLUS(I,J)/ZMAX + 42 CONTINUE + ZMAXC(I) = ZMAX + 30 CONTINUE +*---- +* COST SCALING +*---- + ZMAX = 0.0D0 + DO 45 J=1,N0 + ZMAX = MAX(ZMAX,ABS(COUT(J))) + 45 CONTINUE + DO 50 J=1,N0 + COUT(J) = COUT(J)/ZMAX + 50 CONTINUE +*---- +* PRINT TABLES AFTER SCALING OF COSTS AND CONSTRAINTS +*---- + IF(IMPR.GE.5) THEN + CALL PLNTAB(COUT,APLUS,INPLUS,BGAR,PDG,BINF,BSUP,N0,M0, + > CLNAME//' AFTER SCALING OF COSTS AND CONSTRAINTS') + ENDIF + ALLOCATE(BGAR0(M0)) + DO 52 I=1,M0 + BGAR0(I) = BGAR(I) + 52 CONTINUE +*---- +* INITIAL ESTIMATES +*---- + DELX = SQRT(XDROIT) + EPSS = EPS*DELX +* + DO 55 I=1,N0 + DELTA(I) = DELX/10.0 + RX(I) = 0.0 + 55 CONTINUE + TEMP = DELX/SQRT(REAL(N0))/10 +*---- +* MAP ITERATIONS +*---- + ITER = 0 + 60 ITER = ITER + 1 + CONT = 0.0 + DO 70 I=1,M0+1 + BGAR(I) = BGAR0(I) + 70 CONTINUE +*---- +* CONTROL VARIABLE BOUNDS +*---- + DO 90 I=1,N0 + IF(ITER.EQ.1) THEN +! XOBJ(I) = EPSIR*10.0 + XOBJ(I) = 0.0 + BINF(I) = -TEMP + BSUP(I) = TEMP + ELSE + BINF(I) = -DELTA(I) + BSUP(I) = DELTA(I) + ENDIF +*---- +* LINEARIZATION OF THE QUADRATIC CONSTRAINT +*---- + CONT = CONT + XOBJ(I)**2 + APLUS(M0+1,I) = 2.0*XOBJ(I) + DO 95 J=1,M0 + BGAR(J) = BGAR(J) - APLUS(J,I)*XOBJ(I) + 95 CONTINUE + 90 CONTINUE +* + INPLUS(M0+1) = 1 + BGAR(M0+1) = XDROIT - CONT + M = M0 + 1 +*---- +* REORGANIZE TABLES FOR SIMPLEX +*---- + DO 120 I=1,M + DUMY = 0.0D0 +* + DO 100 J=1,N0 + DUMY = DUMY + APLUS(I,J)*BINF(J) + 100 CONTINUE +* + BGAR(I) = BGAR(I) - DUMY + IF(BGAR(I).GE.0.0) GOTO 120 +* + DO 110 J=1,N0 + APLUS(I,J) = -APLUS(I,J) + 110 CONTINUE +* + BGAR(I) = -BGAR(I) + BGAR0(I)= -BGAR0(I) + BPLUS(I) = -BPLUS(I) + INPLUS(I) = -INPLUS(I) +* + 120 CONTINUE +* + DO 130 J=1,N0 + BSUP(J) = BSUP(J) - BINF(J) + BINF(J) = 0.0 + 130 CONTINUE +*---- +* PRINT SIMPLEX TABLES +*---- + IF(IMPR.GE.5) THEN + CALL PLNTAB(COUT,APLUS,INPLUS,BGAR ,XOBJ ,BINF ,BSUP,N0,M0, + > CLNAME//' AFTER REORGANIZATION FOR SIMPLEX') + ENDIF +*---- +* SOLUTION OF A LINEAR PROGRAMMING PROBLEM USING THE SIMPLEX +*---- + CALL PLSPLX(N0,M,M0+2,1,COUT,APLUS,BGAR,INPLUS,BINF,BSUP,PX, + > DELF,EPSS,IMTHD,IMPR,IERR) +* + DO 140 I=1,N0 + IF(ITER.EQ.1) THEN + PX(I) = PX(I) - TEMP + ELSE + PX(I) = PX(I) - DELTA(I) + ENDIF + 140 CONTINUE +*---- +* SOLUTION OF CURRENT ITERATION +*---- + IF(IMPR.GE.2) THEN + IF(((ITER.GE.1).AND.(IMPR.LE.2)).OR.(IMPR.GE.3)) THEN + WRITE (6,1000) + ENDIF + WRITE (6,2000) ITER,DELF,(PX(I),I=1,N0) + ENDIF +*---- +* DEGENERESCENCE OR EPS TOO SMALL +*---- + IF(IERR.EQ.1) THEN + WRITE(6,3000) ITER + IERR = 3 + RETURN +*---- +* NO SOLUTION IF ITER=1 +*---- + ELSE IF(IERR.EQ.2) THEN + IF(IMPR.GE.1) WRITE(6,4000) ITER + IF(ITER.GE.ITMAX) RETURN + ENDIF +* + ERR = 0.0 + DO 160 I=1,N0 +* + IF((RX(I)*PX(I).LT.0.0).AND.(IERR.EQ.0)) THEN + DELTA(I) = DELTA(I)*0.5 + ENDIF +* + RX(I) = PX(I) + XOBJ(I) = XOBJ(I) + PX(I) + ERR = ERR + PX(I)**2 + 160 CONTINUE +* + ERR = SQRT(ERR) + EPSS = EPS*DELX/10.0 +* + IF(IMPR.GE.1) THEN + WRITE(6,2000) ITER,DELF,(XOBJ(I),I=1,N0) + WRITE(6,2000) ITER,0.0,(DELTA(I),I=1,N0) + ENDIF +* + IF(ERR.LE.EPSS) THEN + IERR = 0 + GOTO 170 + ENDIF +* + IF(ITER.GE.ITMAX) THEN + IERR = 5 + WRITE (6,5000) ITER + RETURN + ENDIF + GO TO 60 +*---- +* RESCALE BACK AND PRINT THE SOLUTION +*---- + 170 DO 175 J=1,N0 + SCAL = SCALE(J) + COUT(J) = COUT(J)*ZMAX/SCAL + XOBJ(J) = XOBJ(J)*SCAL + PDG(J) = XDROIT/SCAL**2 +* + DO 177 I=1,M0 + APLUS(I,J) = APLUS(I,J)/SCAL + 177 CONTINUE + 175 CONTINUE +* + X = 0.0D0 + OBJ = 0.0D0 + DO 180 J=1,N0 + X = X + PDG(J)*XOBJ(J)*XOBJ(J) + OBJ = OBJ + XOBJ(J)*COUT(J) + 180 CONTINUE +* + IF(IMPR.GE.1) THEN + WRITE (6,6000) OBJ,X,(XOBJ(J),J=1,N0) + IF(M0.GT.0) WRITE (6,7000) +* + DO 190 I=1,M0 + XVAL = BPLUS(I) + DO 185 J=1,N0 + XVAL = XVAL - APLUS(I,J)*XOBJ(J)*ZMAXC(I)/SCALE(J) + 185 CONTINUE + WRITE (6,8000) I,XVAL + 190 CONTINUE + ENDIF + DEALLOCATE(ZMAXC,BGAR0) + RETURN +* +1000 FORMAT(/,5X,'ITERATION',8X,'DELF',5X,'CONTROL VARIABLES') +2000 FORMAT(5X,I6,5X,8E12.4,/,(28X,5E12.4)) +3000 FORMAT(5X,I6,5X,'DEGENERESCENCE OR EPS TOO SMALL') +4000 FORMAT(5X,I6,5X,'NO SOLUTION') +5000 FORMAT(5X,I6,5X,'MAXIMUM ITERATION REACHED') +6000 FORMAT(//,5X,'FINAL SOLUTION (MAP1-SIMPLEX) ', + > /,5X,'------------------------', + > /,5X,'OBJECTIVE FUNCTION : ',1P,E12.5, + > /,5X,'QUADRATIC CONSTRAINT : ',1P,E12.5, + > /,5X,'CONTROL VARIABLES : ',/,(10X,10E12.4)) +7000 FORMAT(//,5X,'CONSTRAINT DEVIATIONS : ',/) +8000 FORMAT(2X,I3,'...',2X,1P,D12.4) + END -- cgit v1.2.3