problemfuns.f90

一个基于打靶法的最优控制求解软件 求解过程中采用参数延续算法

Module SpecFuns
  use shootdefs
  implicit none


contains

  Subroutine InitPar(mode, z, lambda)
    implicit none

    integer, intent(in) :: mode
    real(kind=8), intent(in) :: z(n), lambda

    !mode = 0: first call
    !mode = 1: homotopy call

  end Subroutine InitPar


  Subroutine Control(lambda,x,p,u,psi)
    implicit none

    real(kind=8), intent(in) ::  lambda
    real(kind=8), intent(in), dimension(ns) :: x
    real(kind=8), intent(in), dimension(nc) :: p
    real(kind=8), intent(out), dimension(m) :: u
    real(kind=8), intent(out), dimension(dimpsi) :: psi

    !out init
    u = 0d0
    psi = 0d0
    write(0,*) 'Please complete control...'
    stop

  end Subroutine Control

  Subroutine Dynamics(dimphi,lambda,x,p,u,phi)
    implicit none

    integer, intent(in) :: dimphi
    real(kind=8), intent(in) :: lambda
    real(kind=8), intent(in), dimension(ns) :: x
    real(kind=8), intent(in), dimension(nc) :: p
    real(kind=8), intent(in), dimension(m) :: u
    real(kind=8), intent(out), dimension(dimphi) :: phi

    !out init
    phi = 0d0
    write(0,*) 'Please complete dynamics...'
    stop

  end Subroutine Dynamics


  Subroutine FinalConditions(y,s,dsdy)
    implicit none
    
    real(kind=8), dimension(xpdim), intent(in) :: y !!IVP final value
    real(kind=8), dimension(ncf), intent(out) :: s !!Shooting function value   
    real(kind=8), dimension(ncf,xpdim), intent(out) :: dsdy !!Shooting function derivatives
 
    !Local
    integer :: i

    !out init
    s = 0d0
    dsdy = 0d0

    s(1:ncf) = y(fixedT) / yscal(fixedT) - cf(1:ncf)

    do i=1,ncf
       dsdy(i,fixedT(i)) = 1d0 / yscal(fixedT(i))
    end do

  end Subroutine FinalConditions


  Subroutine SwitchFun(lambda,x,p,u,psi)
    implicit none
    real(kind=8), intent(in) :: lambda      
    real(kind=8), intent(in), dimension(ns) :: x         
    real(kind=8), intent(in), dimension(nc) :: p 
    real(kind=8), intent(in), dimension(m) :: u  
    real(kind=8), intent(out), dimension(dimpsi) :: psi  
 
   !out init
    psi = 0d0
    write(0,*) 'Please complete switchfun...'
    stop

  end Subroutine SwitchFun


  Subroutine SwitchGradient(lambda,y,gradpsi)
    implicit none
    real(kind=8), intent(in) :: lambda      
    real(kind=8), intent(in), dimension(ns+nc) :: y         
    real(kind=8), intent(out), dimension(ns+nc) :: gradpsi  
 
    !out init
    gradpsi = 0d0

    gradpsi = gradpsi / yscal(1:ns+nc)
    write(0,*) 'Please complete switchgradient...'
    stop

  end Subroutine SwitchGradient


  Subroutine InitAltcontrol(lambda,alpha,x,p,u)
    implicit none   
    real(kind=8), intent(in) :: lambda      
    real(kind=8), intent(inout) :: alpha  
    real(kind=8), intent(in), dimension(ns) :: x         
    real(kind=8), intent(in), dimension(nc) :: p 
    real(kind=8), intent(inout), dimension(m) :: u  

    !initialize control and alpha
    alpha = (lower(1) + upper(1)) / 2d0


  end Subroutine InitAltcontrol


  Subroutine UpdateAltcontrol(lambda,newalpha,x,p,u)
    implicit none   
    real(kind=8), intent(in) :: lambda      
    real(kind=8), intent(inout) :: newalpha  
    real(kind=8), intent(in), dimension(ns) :: x         
    real(kind=8), intent(in), dimension(nc) :: p 
    real(kind=8), intent(inout), dimension(m) :: u  

    !update control with alpha
    u(1) = newalpha

  end Subroutine UpdateAltcontrol

  Subroutine PreProcess(time,y,x,p)
    implicit none
    real(kind=8), intent(in) :: time
    real(kind=8), intent(in), dimension(ns+nc) :: y
    real(kind=8), intent(out), dimension(ns) :: x
    real(kind=8), intent(out), dimension(nc) :: p

    !out init
    x = 0d0
    p = 0d0

    x(1:ns) = y(1:ns) / yscal(1:ns)
    p(1:nc) = y(ns+1:ns+nc) / yscal(ns+1:ns+nc)

  end Subroutine PreProcess


  Subroutine PostProcess(dim,phi)
    implicit none
    integer, intent(in) :: dim
    real(kind=8), intent(inout), dimension(dim) :: phi

    phi(1:ns) = phi(1:ns) * yscal(1:ns)
    phi(ns+1:ns+nc) = phi(ns+1:ns+nc) * yscal(ns+1:ns+nc)

  end Subroutine PostProcess

SUBROUTINE RHS_DV(dim, time, y, yd, f, fd, lambda, nbdirs)
  IMPLICIT NONE
  INTEGER,INTENT(IN) :: dim
  REAL(KIND=8), DIMENSION(dim),INTENT(OUT) :: f
  REAL(KIND=8),INTENT(IN) :: lambda
  INTEGER :: nbdirs
  REAL(KIND=8) :: fd(nbdirs, dim)
  REAL(KIND=8),INTENT(IN) :: time
  REAL(KIND=8), DIMENSION(dim),INTENT(IN) :: y
  REAL(KIND=8) :: yd(nbdirs, dim)
  write(0,*) 'Please provide the derivatives of the right hand side via automatic differentiation (TAPENADE)'
  stop
end SUBROUTINE RHS_DV

SUBROUTINE SWITCH_DV(dim, ty, tyd, g, gd, nbdirs)
  IMPLICIT NONE
  INTEGER, INTENT(IN) :: dim
  REAL(KIND=8), INTENT(OUT) :: g
  INTEGER :: nbdirs
  REAL(KIND=8), DIMENSION(nbdirs), INTENT(OUT) :: gd
  REAL(KIND=8), DIMENSION(dim+1), INTENT(IN) :: ty
  REAL(KIND=8), DIMENSION(nbdirs, dim+1), INTENT(IN) :: tyd
  INTEGER :: nd

  REAL(KIND=8) :: abs1, abs1d(nbdirs), p(nc), pd(nbdirs, nc)
  REAL(KIND=8) :: time
  REAL(KIND=8) :: x(ns)
  REAL(KIND=8) :: y(dim), yd(nbdirs, dim)
  write(0,*) 'Please provide the derivatives of the switching function via automatic differentiation (TAPENADE)'
  stop
end SUBROUTINE SWITCH_DV


end Module SpecFuns