Skip to content

Commit 1563597

Browse files
jacobwilliamsjacobwilliams
committed
experiment with adding new bounds modes
Fixes #27
1 parent 9f65333 commit 1563597

File tree

2 files changed

+179
-25
lines changed

2 files changed

+179
-25
lines changed

‎src/nlesolver_module.F90‎

Lines changed: 178 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,12 @@ module nlesolver_module
7575
!integer,parameter,public :: NLESOLVER_SPARSITY_CUSTOM_DENSE = 6 !! assume dense (use a user provided dense solver).
7676
! if add will have to update code below where sparsity_modes /= 1 is assumed sparse
7777

78+
! bounds model parameters:
79+
integer,parameter,public :: NLESOLVER_IGNORE_BOUNDS = 0
80+
integer,parameter,public :: NLESOLVER_SCALAR_BOUNDS = 1
81+
integer,parameter,public :: NLESOLVER_VECTOR_BOUNDS = 2
82+
integer,parameter,public :: NLESOLVER_WALL_BOUNDS = 3
83+
7884
!*********************************************************
7985
type,public :: nlesolver_type
8086

@@ -108,12 +114,20 @@ module nlesolver_module
108114
character(len=:),allocatable :: message !! latest status message
109115
integer :: istat = -999 !! latest status message
110116

111-
integer :: bounds_mode = 0 !! how to handle the `x` variable bounds:
117+
integer :: bounds_mode = NLESOLVER_IGNORE_BOUNDS !! how to handle the `x` variable bounds:
112118
!!
113119
!! * 0 = ignore bounds.
114-
!! * 1 = use bounds (if specified) by adjusting the `x` vector
115-
!! at each function evaluation so that each individual `x`
116-
!! component is within its bounds.
120+
!! * 1 = scalar mode : adjust the step vector (`xnew-x`) by moving each individual scalar component
121+
!! of `xnew` to be within the bounds. Note that this can change the direction and magnitude of
122+
!! the search vector.
123+
!! * 2 = vector mode: backtrack the search vector so that no bounds are violated. note that
124+
!! this does not change the direction of the vector, only the magnitude.
125+
!! * 3 = wall mode: adjust the step vector (`xnew-x`) by moving each individual scalar component
126+
!! of `xnew` to be within the bounds. And then holding the ones changed fixed during the line search.
127+
!!
128+
!! Note: if `bounds_mode>0`, then the initial `x` vector is also adjusted
129+
!! to be within the bounds before the first step, if necessary. This is just done
130+
!! by adjusting each scalar component of `x` to be within bounds.
117131
real(wp),dimension(:),allocatable :: xlow !! lower bounds for `x` (size is `n`). only used if `bounds_mode>0`.
118132
real(wp),dimension(:),allocatable :: xupp !! upper bounds for `x` (size is `n`). only used if `bounds_mode>0`.
119133

@@ -177,6 +191,8 @@ module nlesolver_module
177191

178192
procedure :: set_status
179193
procedure :: adjust_x_for_bounds
194+
procedure :: adjust_search_direction
195+
procedure :: compute_next_step
180196

181197
end type nlesolver_type
182198
!*********************************************************
@@ -354,6 +370,7 @@ end function real2str
354370
! * -14 -- Error: must specify grad_sparse, irow, and icol for sparsity_mode > 1
355371
! * -15 -- Error: irow and icol must be the same length
356372
! * -16 -- Error: xlow > xupp
373+
! * -17 -- Error adjusting line search direction for bounds
357374
! * -999 -- Error: class has not been initialized
358375
! * 0 -- Class successfully initialized in [[nlesolver_type:initialize]]
359376
! * 1 -- Required accuracy achieved
@@ -475,9 +492,9 @@ subroutine initialize_nlesolver_variables(me,&
475492
integer,intent(in),optional :: bounds_mode !! how to handle the `x` variable bounds:
476493
!!
477494
!! * 0 = ignore bounds
478-
!! * 1 = use bounds (if specified) by adjusting the `x` vector
479-
!! at each step so that each individual `x` component is within
480-
!! the bounds
495+
!! * 1 = scalar mode
496+
!! * 2 = vector mode
497+
!! * 3 = wall mode
481498
real(wp),dimension(n),intent(in),optional :: xlow !! lower bounds for `x` (size is `n`). only used if `bounds_mode>0` and
482499
!! both `xlow` and `xupp` are specified.
483500
real(wp),dimension(n),intent(in),optional :: xupp !! upper bounds for `x` (size is `n`). only used if `bounds_mode>0` and
@@ -508,7 +525,7 @@ subroutine initialize_nlesolver_variables(me,&
508525
me%xupp = xupp
509526
me%xlow = xlow
510527
else
511-
me%bounds_mode = 0 ! default
528+
me%bounds_mode = NLESOLVER_IGNORE_BOUNDS ! default
512529
end if
513530

514531
if (present(step_mode)) then
@@ -985,6 +1002,7 @@ end subroutine nlesolver_solver
9851002
!*****************************************************************************************
9861003
!>
9871004
! if necessary, adjust the `x` vector to be within the bounds.
1005+
! used for the initial guess.
9881006

9891007
subroutine adjust_x_for_bounds(me,x)
9901008

@@ -995,22 +1013,138 @@ subroutine adjust_x_for_bounds(me,x)
9951013

9961014
integer :: i !! counter
9971015

998-
if (me%bounds_mode==1) then
1016+
if (me%bounds_mode/=NLESOLVER_IGNORE_BOUNDS) then
1017+
! for all bounds modes, adjust the initial guess to be within the bounds
9991018
! x = min(max(x,me%xlow),me%xupp)
10001019
do i = 1, me%n
10011020
if (x(i)<me%xlow(i)) then
10021021
x(i) = me%xlow(i)
1003-
if (me%verbose) write(me%iunit, '(A)') 'x('//int2str(i)//') < xlow(i) : adjusting to lower bound'
1022+
if (me%verbose) write(me%iunit, '(A)') 'Initial x('//int2str(i)//') < xlow(i) : adjusting to lower bound'
10041023
else if (x(i)>me%xupp(i)) then
10051024
x(i) = me%xupp(i)
1006-
if (me%verbose) write(me%iunit, '(A)') 'x('//int2str(i)//') > xupp(i) : adjusting to upper bound'
1025+
if (me%verbose) write(me%iunit, '(A)') 'Initial x('//int2str(i)//') > xupp(i) : adjusting to upper bound'
10071026
end if
10081027
end do
10091028
end if
10101029

10111030
end subroutine adjust_x_for_bounds
10121031
!*****************************************************************************************
10131032

1033+
!*****************************************************************************************
1034+
!>
1035+
! if necessary, adjust the search direction vector so that bounds are not violated.
1036+
!
1037+
!### Reference
1038+
! * https://openmdao.org/newdocs/versions/latest/features/building_blocks/solvers/bounds_enforce.html
1039+
1040+
! TODO: for the 'wall' mode, flag the ones that were violated, and use
1041+
! that to set the alpha=0 for those vars during the line search.
1042+
! --> where (adjusted) xnew = x
1043+
1044+
subroutine adjust_search_direction(me,x,p,pnew,modified)
1045+
1046+
implicit none
1047+
1048+
class(nlesolver_type),intent(inout) :: me
1049+
real(wp),dimension(me%n),intent(in) :: x !! initial `x`
1050+
real(wp),dimension(me%n),intent(in) :: p !! search direction `p = xnew - x`
1051+
real(wp),dimension(me%n),intent(out) :: pnew !! new search direction
1052+
logical,dimension(me%n),intent(out) :: modified !! indicates the elements of p that were modified
1053+
1054+
integer :: i !! counter
1055+
real(wp),dimension(:),allocatable :: xnew !! `x + pnew`
1056+
logical :: search_direction_modifed !! if `p` has been modified
1057+
real(wp) :: t !! indep var in the step equation: `xnew = x + t*p`
1058+
1059+
if (me%bounds_mode==NLESOLVER_IGNORE_BOUNDS) then ! no change
1060+
modified = .false.
1061+
pnew = p
1062+
return
1063+
end if
1064+
1065+
allocate(xnew(me%n))
1066+
xnew = x + p
1067+
search_direction_modifed = .false.
1068+
t = 1.0_wp ! for mode 2, start with full newton step
1069+
! note: have to check each variable and choose the smallest t.
1070+
modified = .false. ! initialize
1071+
1072+
do i = 1, me%n
1073+
if (xnew(i)<me%xlow(i)) then
1074+
search_direction_modifed = .true.
1075+
modified(i) = .true.
1076+
if (me%verbose) write(me%iunit, '(A)') 'x('//int2str(i)//') < xlow(i) : adjusting to lower bound'
1077+
select case (me%bounds_mode)
1078+
case(NLESOLVER_SCALAR_BOUNDS,NLESOLVER_WALL_BOUNDS)
1079+
xnew(i) = me%xlow(i)
1080+
case(NLESOLVER_VECTOR_BOUNDS)
1081+
t = min(t,(me%xlow(i)-x(i))/p(i))
1082+
end select
1083+
else if (xnew(i)>me%xupp(i)) then
1084+
search_direction_modifed = .true.
1085+
modified(i) = .true.
1086+
if (me%verbose) write(me%iunit, '(A)') 'x('//int2str(i)//') > xupp(i) : adjusting to upper bound'
1087+
select case (me%bounds_mode)
1088+
case(NLESOLVER_SCALAR_BOUNDS,NLESOLVER_WALL_BOUNDS)
1089+
xnew(i) = me%xupp(i)
1090+
case(NLESOLVER_VECTOR_BOUNDS)
1091+
t = min(t,(me%xupp(i)-x(i))/p(i))
1092+
end select
1093+
end if
1094+
end do
1095+
1096+
if (search_direction_modifed) then
1097+
! adjust the search direction:
1098+
select case (me%bounds_mode)
1099+
case (NLESOLVER_SCALAR_BOUNDS,NLESOLVER_WALL_BOUNDS)
1100+
pnew = xnew - x ! here we have changed the search direction vector
1101+
if (all(pnew==0.0_wp)) call me%set_status(istat = -17, string = 'Error adjusting line search direction for bounds')
1102+
case (NLESOLVER_VECTOR_BOUNDS)
1103+
! here was are staying on the original search direction vector, just walking back
1104+
if (t <= 0.0_wp) then ! something wrong
1105+
call me%set_status(istat = -17, string = 'Error adjusting line search direction for bounds')
1106+
pnew = p
1107+
else
1108+
pnew = p / t
1109+
end if
1110+
end select
1111+
if (me%verbose) write(me%iunit, '(A)') 'Search direction modified to be within bounds'
1112+
else
1113+
pnew = p
1114+
end if
1115+
1116+
end subroutine adjust_search_direction
1117+
!*****************************************************************************************
1118+
1119+
!*****************************************************************************************
1120+
!>
1121+
! Compute the next step.
1122+
1123+
subroutine compute_next_step(me, xold, search_direction, alpha, modified, xnew)
1124+
1125+
class(nlesolver_type),intent(inout) :: me
1126+
real(wp),dimension(me%n),intent(in) :: xold !! initial `x`
1127+
real(wp),dimension(me%n),intent(in) :: search_direction !! search direction vector
1128+
real(wp),intent(in) :: alpha !! step length
1129+
logical,dimension(me%n),intent(in) :: modified !! indicates the elements of `p` that were
1130+
!! modified because they violated the bounds
1131+
real(wp),dimension(me%n),intent(out) :: xnew !! final `x`
1132+
1133+
if (me%bounds_mode == NLESOLVER_WALL_BOUNDS) then
1134+
! for the 'wall' mode, the modified variables are held fixed during the step
1135+
where (modified)
1136+
xnew = xold
1137+
else where
1138+
xnew = xold + search_direction * me%alpha
1139+
end where
1140+
else
1141+
! all other modes just use the computed search direction
1142+
xnew = xold + search_direction * alpha
1143+
end if
1144+
1145+
end subroutine compute_next_step
1146+
!*****************************************************************************************
1147+
10141148
!*****************************************************************************************
10151149
!>
10161150
! Destructor
@@ -1129,8 +1263,14 @@ subroutine simple_step(me,xold,p,x,f,fvec,fjac,fjac_sparse)
11291263
real(wp),dimension(:,:),intent(in),optional :: fjac !! jacobian matrix [dense]
11301264
real(wp),dimension(:),intent(in),optional :: fjac_sparse !! jacobian matrix [sparse]
11311265

1132-
x = xold + p * me%alpha
1133-
call me%adjust_x_for_bounds(x)
1266+
real(wp),dimension(:),allocatable :: search_direction
1267+
logical,dimension(:),allocatable :: modified !! indicates the elements of p that were modified
1268+
1269+
allocate(search_direction(me%n))
1270+
allocate(modified(me%n))
1271+
1272+
call me%adjust_search_direction(xold,p,search_direction,modified)
1273+
call me%compute_next_step(xold, search_direction, me%alpha, modified, x)
11341274

11351275
!evaluate the function at the new point:
11361276
call me%func(x,fvec)
@@ -1170,11 +1310,18 @@ subroutine backtracking_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
11701310
real(wp),dimension(:),allocatable :: gradf !! line search objective function gradient vector
11711311
real(wp),dimension(:),allocatable :: xtmp !! `x` value for linesearch
11721312
real(wp),dimension(:),allocatable :: fvectmp !! `fvec` value for linesearch
1313+
real(wp),dimension(:),allocatable :: search_direction !! search direction to use (can be modified from `p` if bounds are violated)
1314+
logical,dimension(:),allocatable :: modified !! indicates the elements of p that were modified
11731315

11741316
! allocate arrays:
11751317
allocate(gradf(me%n))
11761318
allocate(xtmp(me%n))
11771319
allocate(fvectmp(me%m))
1320+
allocate(search_direction(me%n))
1321+
allocate(modified(me%n))
1322+
1323+
! set the search direction:
1324+
call me%adjust_search_direction(xold,p,search_direction,modified)
11781325

11791326
! compute the gradient of the function to be minimized
11801327
! (which in this case is 1/2 the norm of fvec). Use the chain
@@ -1190,7 +1337,7 @@ subroutine backtracking_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
11901337
gradf(i) = dot_product(fvec,pack(fjac_sparse,mask=me%icol==i))
11911338
end do
11921339
end if
1193-
slope = dot_product(p, gradf)
1340+
slope = dot_product(search_direction, gradf)
11941341
t = -me%c * slope
11951342

11961343
if (me%verbose) then
@@ -1199,12 +1346,12 @@ subroutine backtracking_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
11991346
end if
12001347

12011348
! perform the line search:
1202-
min_alpha_reached = .false.
1203-
alpha = me%alpha_max ! start with the largest step
1349+
1350+
min_alpha_reached = .false. ! initialize
1351+
alpha = me%alpha_max ! start with the largest step
12041352
do
12051353

1206-
xtmp = xold + p * alpha
1207-
call me%adjust_x_for_bounds(xtmp)
1354+
call me%compute_next_step(xold, search_direction, alpha, modified, xtmp)
12081355
call me%func(xtmp,fvectmp)
12091356
ftmp = norm2(fvectmp)
12101357

@@ -1264,16 +1411,20 @@ subroutine exact_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
12641411

12651412
real(wp),dimension(:),allocatable :: xnew !! used in [[func_for_fmin]]
12661413
real(wp) :: alpha_min
1414+
real(wp),dimension(:),allocatable :: search_direction !! search direction to use (may be modified from `p` if bounds are violated)
1415+
logical,dimension(:),allocatable :: modified !! indicates the elements of p that were modified
12671416

12681417
allocate(xnew(me%n))
1418+
allocate(search_direction(me%n))
1419+
allocate(modified(me%n))
12691420

12701421
! find the minimum value of f in the range of alphas:
12711422
alpha_min = fmin(func_for_fmin,me%alpha_min,me%alpha_max,me%fmin_tol)
12721423

12731424
if (me%verbose) write(me%iunit,'(1P,*(A,1X,E16.6))') ' alpha_min = ', alpha_min
12741425

1275-
x = xold + p * alpha_min
1276-
call me%adjust_x_for_bounds(x)
1426+
call me%adjust_search_direction(xold,p,search_direction,modified)
1427+
call me%compute_next_step(xold, search_direction, alpha_min, modified, x)
12771428
if (all(x==xnew)) then
12781429
! already computed in the func
12791430
else
@@ -1288,8 +1439,7 @@ real(wp) function func_for_fmin(alpha)
12881439
implicit none
12891440
real(wp),intent(in) :: alpha !! indep variable
12901441

1291-
xnew = xold + p * alpha
1292-
call me%adjust_x_for_bounds(xnew)
1442+
call me%compute_next_step(xold, search_direction, alpha, modified, xnew)
12931443
call me%func(xnew,fvec)
12941444
func_for_fmin = norm2(fvec) ! return result
12951445

@@ -1326,6 +1476,8 @@ subroutine fixed_point_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
13261476
real(wp) :: f_tmp !! temp `f`
13271477
real(wp) :: step_size !! step size for `alpha`
13281478
integer :: n !! number of steps to divide the interval
1479+
real(wp),dimension(:),allocatable :: search_direction !! search direction to use (may be modified from `p` if bounds are violated)
1480+
logical,dimension(:),allocatable :: modified !! indicates the elements of p that were modified
13291481

13301482
! 1 o-----------o
13311483
! 2 o-----o-----o
@@ -1337,6 +1489,8 @@ subroutine fixed_point_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
13371489
allocate(alphas_to_try(n_points))
13381490
allocate(x_tmp(me%n))
13391491
allocate(fvec_tmp(me%m))
1492+
allocate(search_direction(me%n))
1493+
allocate(modified(me%n))
13401494

13411495
step_size = (me%alpha_max - me%alpha_min) / real(n,wp)
13421496

@@ -1349,10 +1503,10 @@ subroutine fixed_point_linesearch(me,xold,p,x,f,fvec,fjac,fjac_sparse)
13491503

13501504
! now compute the functions at these alphas:
13511505
f = big
1506+
call me%adjust_search_direction(xold,p,search_direction,modified)
13521507
do i = 1, n_points
13531508

1354-
x_tmp = xold + p * alphas_to_try(i)
1355-
call me%adjust_x_for_bounds(x_tmp)
1509+
call me%compute_next_step(xold, search_direction, alphas_to_try(i), modified, x_tmp)
13561510

13571511
! evaluate the function at tthis point:
13581512
call me%func(x_tmp,fvec_tmp)

‎test/nlesolver_test_1.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -99,7 +99,7 @@ program nlesolver_test_1
9999
n_intervals = n_intervals, &
100100
fmin_tol = fmin_tol, &
101101
verbose = verbose, &
102-
bounds_mode = 1, &
102+
bounds_mode = NLESOLVER_SCALAR_BOUNDS, &
103103
xlow = [0.0_wp, -5.0_wp], &
104104
xupp = [1.0_wp, 0.0_wp])
105105
call solver%status(istat, message)

0 commit comments

Comments
 (0)