multiplier.f90

求解约束优化问题增广拉格朗日乘子法源程序

module Multipliter_Variable
parameter(n_dimension=100)
implicit doubleprecision (a-h,o-z)

double precision::u0(n_dimension),v0(n_dimension)
double precision::a0,b0,c0,error,ph

integer::l,m,lh,lf,lj
logical::show

end module
subroutine Multipliter_distribution(n,xx,ii)
use Multipliter_Variable
implicit doubleprecision (a-h,o-z)
double precision::xx(n_dimension),hx(n_dimension),qx(n_dimension)

lj=ii
call Initialization()
call Phr(n,xx,fx,hx,qx)
end subroutine

subroutine Multipliter(n,xx,ii)
use Multipliter_Variable
implicit doubleprecision (a-h,o-z)
double precision::xx(n_dimension),hx(n_dimension),qx(n_dimension)

lj=ii
call Initialization()
call Flow_mass2Pressure_drop(n,xx,fx,hx,qx) 
end subroutine


subroutine Initialization()
use Multipliter_Variable
implicit doubleprecision (a-h,o-z)
double precision::xx(n_dimension),hx(n_dimension),qx(n_dimension)

l=1
m=0
a0=4
b0=0.5
c0=8
error=1.0e-2
show=.TRUE.
do i=1,l
    u0(i)=1.0
end do
do i=1,m
    v0(i)=1.0
end do
!allocate(xx(n),hx(n),qx(n),u0(n),v0(n)) 

end subroutine


subroutine Phr(n,xx,fx,hx,qx)
use Multipliter_Variable
implicit doubleprecision (a-h,o-z)
double precision::xx(n_dimension),hx(n_dimension),qx(n_dimension)
double precision::u1(n_dimension)

!allocate(xx(n),hx(n),qx(n),u1(n)) 
lh=1
lf=1
call ObjFun(n,xx,fx,hx,qx)
ph0=0.0
do i=1,l
    ph0=ph0+hx(i)*hx(i)
end do
ph0=sqrt(ph0)
do while(.true.)
    call Hook(n,xx,fx,hx,qx)
	if(abs(ph)<error) exit
    if(show) then
!       write(*,(5X,'lh=',I2,5X,'lf=',I6,5X,'ph=',F20.10)) lh,lf,ph
!		write(*,(5X,'f(xx)=',F20.10)) fx
!		write(*,(6X,'xx(',i2,')=',F20.10)) (i,xx(i),i=1,n)
!		write(*,(6X,'h(',i2,')=',F20.10)) (i,h(i),i=1,n)
!		write(*,(6X,'q(',i2,')=',F20.10)) (i,q(i),i=1,n)
        write(*,*) lh,lf,ph
		write(*,*) fx
		write(*,*) (i,xx(i),i=1,n)
		write(*,*) (i,hx(i),i=1,n)
		write(*,*) (i,qx(i),i=1,n)
	end if
	ph1=ph/ph0
	ph0=ph
	if((ph1-b0)<0) then
	    c0=1.0*c0
	else
	    c0=a0*c0
	end if 
	do i=1,l
        u1(i)=u0(i)-c0*hx(i)
		u0(i)=u1(i)
	end do
	do i=1,m
        ab=v0(i)-c0*qx(i)
		v0(i)=(abs(ab)+ab)/2.0
	end do
    lh=lh+1
end do
return
end subroutine


subroutine Prx(n,xx,fx,hx,qx,px)
use Multipliter_Variable
implicit doubleprecision (a-h,o-z)
double precision::xx(n_dimension),hx(n_dimension),qx(n_dimension)

!allocate(xx(n),hx(n),qx(n)) 
call ObjFun(n,xx,fx,hx,qx)
uh=0.0
hh=0.0
do i=1,l
    uh=uh+u0(i)*hx(i)
	hh=hh+hx(i)*hx(i)
end do
ph=sqrt(hh)
pq=0.0
do i=1,m
    aa=v0(i)-c0*qx(i)
    pq=pq+((abs(aa)+aa)/2.0)**2-v0(i)*v0(i)
end do
px=fx-uh+(0.5*c0)*hh+(0.5/c0)*pq
return
end subroutine


subroutine Hook(n,xx,fx,hx,qx)
use Multipliter_Variable
implicit doubleprecision (a-h,o-z)
double precision::xx(n_dimension),hx(n_dimension),qx(n_dimension)
double precision::y(n_dimension),z(n_dimension)

!allocate(xx(n),hx(n),qx(n),y(n),z(n)) 
h=0.1
t0=1.0
do while(.true.)
    do i=1,n
        y(i)=xx(i)
    end do
	do while(.true.)
        do i=1,n
            call Prx(n,xx,fx,hx,qx,p1)
			call Prx(n,y,fx,hx,qx,p2)
			y(i)=y(i)+h
			call Prx(n,y,fx,hx,qx,p3)
			if(p3>=p2) then
                y(i)=y(i)-2.0*h
				call Prx(n,y,fx,hx,qx,p4)
				if(p4>=p2) then
                    y(i)=y(i)+h
				end if
			end if
		end do
		call Prx(n,y,fx,hx,qx,p5)
		if(p5>=p1) exit
		do i=1,n
		    z(i)=y(i)
			y(i)=z(i)+t0*(z(i)-xx(i))
			xx(i)=z(i)
		end do
	end do
	if(h<error) exit
	h=0.5*h
end do
end subroutine


subroutine ObjFun(n,xx,fx,hx,qx)
use Multipliter_Variable
implicit doubleprecision (a-h,o-z)
double precision::xx(n_dimension),hx(n_dimension),qx(n_dimension),mass

!please inter you function        

end subroutine