Subroutine optsensitivity
633 :
634 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
635 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
636 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
637 : ! ADDED BY STAN TO ASSESS SENSITIVITY OF OBJECTIVE FUNCTION TO EACH OPTIMISED
638 : ! PARAMETER:
639 : !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
640 :
641 : subroutine optsensitivity ()
642 : use vom_sce_mod
643 : implicit none
644 :
645 : INTEGER :: ii, jj, kk, failed10, pos
646 : REAL*8 :: distmin, distmax
647 : REAL*8 :: oldpar, newpar
648 : REAL*8 :: ofvec2, ofchange, parchange
649 : CHARACTER(300) :: writeformat
650 : CHARACTER(len=135) :: msg
651 : REAL*8, ALLOCATABLE :: tmp_8(:)
652 :
653 : shufflevar2(:) = shufflevar(:,1)
654 : ofvec2 = ofvec(1)
655 : failed10 = 0
656 :
657 : write(kfile_progress,*) 'SENSITIVITY ANALYSIS'
658 : write(kfile_progress,*) 'changes in % of feasible range)'
659 : dataarray(1,1:nopt) = shufflevar2(optid(:))
660 : dataarray(1,nopt+1) = ofvec2
661 : evolution = 'test'
662 : pos = 0
663 :
664 : do ii = 1, nopt
665 :
666 : oldpar = shufflevar2(optid(ii))
667 :
668 : ! * TO MAKE SURE THAT PERTURBATIONS DO NOT EXCEED THE FEASIBLE RANGE
669 : ! * AND THE PARAMETER VALUES THEMSELVES
670 :
671 : distmin = oldpar - parmin(optid(ii))
672 : distmax = parmax(optid(ii)) - oldpar
673 : write(kfile_progress,*) " "
674 : writeformat = '("change of var: ",A9, "(",E9.3,")")'
675 : write(msg,writeformat) parname(optid(ii)), oldpar
676 : write(kfile_progress,*) TRIM(msg)
677 : writeformat = '(f7.3,"% (",e14.7,")",": change of OF by ",'
678 : writeformat(44:68) = 'e10.3,"%"," (",e10.3,")")'
679 :
680 : jj = 0 ! counts from 0 to 3 and than from 3 to 0
681 : kk = 1 ! counter: first +1 after 3 steps it becomes -1
682 :
683 : do while (jj .ge. 0)
684 :
685 : if (kk .eq. 1) then
686 : newpar = oldpar - distmin * 10.d0 ** (-jj)
687 : else
688 : newpar = oldpar + distmax * 10.d0 ** (-jj)
689 : endif
690 : nrun = nrun + 1
691 : shufflevar2(optid(ii)) = newpar
692 :
693 : call transpmodel(shufflevar2(:), SIZE(shufflevar2(:)), ofvec2, 1)
694 :
695 : if (ofvec2 .gt. bestobj) then
696 : bestobj = ofvec2
697 : call write_lastbest(shufflevar2(:), vom_npar, bestobj, 0)
698 : endif
699 :
700 : ! * use temporary variable to prevent warning in ifort
701 : allocate(tmp_8(nopt))
702 : tmp_8(:) = shufflevar2(optid(:))
703 : write(kfile_sceout,outformat) tmp_8(:), ofvec2
704 : deallocate(tmp_8)
705 :
706 : if (kk .eq. 1) then
707 : dataarray((ii-1)*8+jj+2,1:nopt) = shufflevar2(optid(:))
708 : dataarray((ii-1)*8+jj+2,nopt+1) = ofvec2
709 : else
710 : dataarray((ii-1)*8+jj+6,1:nopt) = shufflevar2(optid(:))
711 : dataarray((ii-1)*8+jj+6,nopt+1) = ofvec2
712 : endif
713 : ofchange = (ofvec2 - dataarray(1,nopt + 1)) &
714 : & / ABS(dataarray(1,nopt + 1)) * 100.d0
715 : parchange = (newpar - oldpar) / (parmax(optid(ii)) &
716 : & - parmin(optid(ii))) * 100.d0 ! the change of the parameter in % of feasible range
717 : write(msg,writeformat) parchange, newpar, ofchange, ofvec2
718 : write(kfile_progress,*) TRIM(msg)
719 : flush(kfile_progress)
720 :
721 : if (ofchange .gt. 1.0d-10) then
722 : shufflevar(:,sopt-pos) = shufflevar2(:)
723 : ofvec(sopt-pos) = ofvec2
724 : pos = pos + 1
725 : if (ABS(parchange) .gt. i_resolution) then
726 : failed10 = failed10 + 1
727 : endif
728 : endif
729 :
730 : jj = jj + kk
731 : if (jj .eq. 4) then
732 : jj = 3
733 : kk = -1
734 : endif
735 : enddo
736 :
737 : shufflevar2(optid(ii)) = oldpar
738 :
739 : enddo
740 :
741 : ! * ASSESS SECOND CONVERGENCE CRITERIUM: NO PARAMETER CHANGE BY MORE THAN
742 : ! * 10% LEADS TO AN INCREASE IN OBJECTIVE FUNCTION
743 :
744 : if (failed10 .gt. 0) then
745 : writeformat = '(I2," parameter(s) more than ",F6.3,"% out of optimum.")'
746 : write(msg,writeformat) failed10, i_resolution
747 : write(kfile_progress,*) TRIM(msg)
748 : write(kfile_progress,*) "Optimisation continued..."
749 : else
750 : write(kfile_progress,*) " "
751 : write(kfile_progress,*) " "
752 : write(kfile_progress,*) "Second convergence criterion satisfied:"
753 : writeformat = '("no parameter shift by more than ",F6.3,"% of max distance leads to")'
754 : write(msg,writeformat) i_resolution
755 : write(kfile_progress,*) TRIM(msg)
756 : write(kfile_progress,*) "an increase in the objective function."
757 : success = 1
758 : endif
759 :
760 : ! * IF CRITERIUM NOT SATISFIED, APPEND NEWLY CREATED DATA TO THE END OF
761 : ! * shufflevar AND ofvec ARRAYS AND RETURN TO MAIN LOOP TO RE-RUN SCE-LOOP
762 : ! * --> deleted
763 :
764 : return
765 : end subroutine optsensitivity