On this page, I will attempt to guide you through the various routines in order to create a 'picture' for you to understand what's going on. First, the code required and then a sample of output summary.
1. Introduce the Find statement that is used to tweak variables until ones objective(s) are met. The variables may have bounds and constraints that restrick where the selected solver may search for an Optimal solution.
! Find statement is key to tweak variables AND obtain an Optimal
! (i.e., Maximum / Minimum) solutions.
FIND gain, Yzero, Preal, Pimag; ! variable list of those to tweak!
in fitboth; by Jove; ! Jove is the solver for this problem
with lower h8low; and uppers h8hi; ! lower & upper bounds
Holding circleUp, PoleWide, HoldZero, ZeroWide; ! constraints
TO MINIMIZE errsum ! objective Minimize errSum variable!
2. This is a sample of code necessary to solve this problem, but hope it gives you an idea what is going on. Variables to tweak (e.g. gain) MUST always be on the rightside of an equal sign; so the solver can change a value when it wants too. Follow gain variable throughout the code for seeing what is required. Plus, see gain and other tweaked variables in output section.
model fitboth
include 'freq68.inc'
C ... Zero limits ... >= 0
do 20 j = 1, nYzero
20 holdZero(j) = (Yzero-Zbottom)
if( nZpairs + nZquads .gt. 0) then
do 22 j = nYzero+1, nYzero + nZpairs + nZquads
ZeroWide(j) = (fnorm - Zimag(j))
22 holdZero(j) = (Zreal(j)-Zbottom) * (Xlmax - Zimag(j))
endif
C ... Pole limits ... >= 0
if( nYpole .gt. 0) then
do 30 j = 1, nYpole
PoleWide(j) = (fnorm - Ypole)
30 circleUp(j) = circle( Ypole, 0, radius)
endif
if( nPpairs .gt. 0) then
do 32 j = 1, nPpairs
PoleWide(j+nYpole) = (fnorm - Pimag(j))
32 circleUp(j+nYpole) = circle( Preal(j), Pimag(j), radius)
endif
call transfer ! calc. 'errMag' error
call Gdelay ! calc. 'rip8sq' error
errMag = 1D10 * errMag * wt8magn
errsum = errMag + rip8sq * 10 ! errSum, is it decreasing in amplitude?
end
Fmodel circle( Vreal, Vimag, radius) ! constrain Vreal & Vimag to within radius?
circle = radius - sqrt(Vreal**2 + Vimag**2)
end
model transfer
include 'freq68.inc'
errMag = 0.
do 50 ii = 1, npoints ! --- CALCULATE TRANSFER FUNCTION ----
F = freq( ii): Hw = Hs( F)
error( ii) = Hw * y8in( ii) - y8out( ii) ! Absolute Error
error( ii) = error( ii) * (y8out( ii)**2) ! Relative
errMag = errMag + error( ii)**2
50 continue
end
Fmodel Hs( F)
include 'freq68.inc'
real*8 num
F2 = F * F: num = 1 : den = 1
if( nYpole + nYpoleHo .gt. 0) den = factor( F2, -Ypole, 0.)
if( nPholds .gt. 0) then
do 20 ij = 1, nPholds
den = den * factor( F2, -PrealH( ij), PimagH( ij))
20 continue
endif
if( nPpairs .gt. 0) then
do 24 ij = 1, nPpairs
den = den * factor( F2, -Preal( ij), Pimag( ij))
24 continue
endif
if( nYzero .gt. 0) then
num = num * factor( F2, -Yzero, 0.)
if( doublet .gt. 0) num = num**2
if( lossless .gt. 0) num = num**2
endif
if( nXzeros .gt. 0) then
num = 100 * num
do 40 ij = 1, nXzeros
num = num * factor( F2, 0., Xzeros(ij))
40 continue
endif
if( nZpairs .gt. 0) then
do 60 ij = 1, nZpairs
num = num * factor( F2, -Zreal(ij), Zimag( ij))
60 continue
endif
n = nZpairs
if( nZquads .gt. 0) then
do 80 ij = 1, nZquads
n = n + 1
num = num * factor( F2, Zreal( n), Zimag( n))
num = num * factor( F2, -Zreal( n), Zimag( n))
80 continue
endif
q = gain * sqrt( abs(num / den))
if( q .gt. 1.D20) q = 1.D20
Hs = q
end
Fmodel factor( F2, sigma, omega)
r2 = sigma**2
if( omega .eq. 0.D0) then
factor = 1. : if( sigma .eq. 0.D0) return ! not sure on value
factor = (r2 - F2) / 10 ! '10' normalizing factor
return
endif
o2 = omega**2 : sum = (r2 + o2 - F2) / 10
temp = sum*sum + 4*F2*r2/100 : factor = 0
if( temp .eq. 0.D0) return
temp = sqrt( temp)
factor = temp / (r2 + o2) ! this R2+O2 is 4 normalizing pole values
end
ooo ! and adjusts Gain value for system.
3. The following block shows a sample of output from the JOVE solver. The first column lists the variables being tweaked; some variables are arrays. The next column, INITIAL, shows initial values for variables that gets things rolling. The next columns are the values for that iteration. This problem required 13 iterations to solve for an optimal solution. This is the most complicated problem I solved todate. Solution time was less than 4 minute on a 5-year old PC.
--- JOVE SUMMARY, INVOKED AT FIT[98] FOR MODEL FITBOTH ----
CONVERGENCE CONDITION AFTER 13 ITERATIONS
OBJECTIVE CRITERION SATISFIED
ALL SPECIFIED CRITERIA SATISFIED
LOOP NUMBER ......... [INITIAL] 1 2
UNKNOWNS
gain 2.500000E-05 1.821383E-01 1.822404E-01
Yzero 5.000000E-01 9.857444E-01 9.861072E-01
Preal ( 1) 2.000000E-01 7.472368E-01 7.469582E-01
Preal ( 2) 2.000000E-01 5.662083E-01 5.651061E-01
Preal ( 3) 2.000000E-01 3.579208E-02 3.750613E-02
Preal ( 4) 2.000000E-01 2.776521E-01 2.782288E-01
PIMAG ( 1) 5.000000E-03 6.431970E-02 6.747955E-02
PIMAG ( 2) 2.800000E-01 4.918416E-01 4.931076E-01
PIMAG ( 3) 5.000000E-01 7.491454E-01 7.490616E-01
PIMAG ( 4) 5.000000E-01 5.795703E-01 5.808899E-01
OBJECTIVE
ERRSUM 1.617663E+07 1.143457E+03 1.101378E+03
INEQUALITY CONSTRAINTS
CIRCLEUP( 1) 5.499375E-01 7.444890E-08 1.745802E-10
CIRCLEUP( 2) 4.059070E-01 1.840798E-08 4.869916E-11
CIRCLEUP( 3) 2.114835E-01 2.105579E-08 5.069178E-11
CIRCLEUP( 4) 2.114835E-01 1.073551E-01 1.059159E-01
POLEWIDE( 1) 1.999500E+01 1.993568E+01 1.993252E+01
POLEWIDE( 2) 1.972000E+01 1.950816E+01 1.950689E+01
POLEWIDE( 3) 1.950000E+01 1.925085E+01 1.925094E+01
POLEWIDE( 4) 1.950000E+01 1.942043E+01 1.941911E+01
ooo
LOOP NUMBER ......... [INITIAL] 13
UNKNOWNS
gain 2.500000E-05 1.822269E-01
Yzero 5.000000E-01 9.861653E-01
Preal ( 1) 2.000000E-01 7.469539E-01
Preal ( 2) 2.000000E-01 5.650895E-01
Preal ( 3) 2.000000E-01 3.753097E-02
Preal ( 4) 2.000000E-01 2.782460E-01
PIMAG ( 1) 5.000000E-03 6.752736E-02
PIMAG ( 2) 2.800000E-01 4.931266E-01
PIMAG ( 3) 5.000000E-01 7.490603E-01
PIMAG ( 4) 5.000000E-01 5.809214E-01
OBJECTIVE
ERRSUM 1.617663E+07 1.100697E+03 ! always decreasing?
INEQUALITY CONSTRAINTS
CIRCLEUP( 1) 5.499375E-01 9.658940E-15
CIRCLEUP( 2) 4.059070E-01 1.332268E-14
CIRCLEUP( 3) 2.114835E-01 4.606050E-08
CIRCLEUP( 4) 2.114835E-01 1.058801E-01
POLEWIDE( 1) 1.999500E+01 1.993247E+01
POLEWIDE( 2) 1.972000E+01 1.950687E+01
POLEWIDE( 3) 1.950000E+01 1.925094E+01
POLEWIDE( 4) 1.950000E+01 1.941908E+01
---END OF LOOP SUMMARY
Are you ready to give it a try? First step is to download our free Calculus (Level) Compiler (i.e., FC-Compiler), install it, and run some Demo (example) files. Next, copy another (Demo) file and save it under another name in your 'user' folder. Edit this file with your equations and (math) models, then save it. Next, try executing it and see how it goes. Enjoy!