* ------------------------------------------------------------------ * File spring.f for SNOPT * * This is a main program to generate an optimal control problem * of arbitrary size and solve it by calling SNOPT as a subroutine. * * The problem size depends on a parameter T. There are * 2T constraints and 3T + 2 variables, as well as bounds * on the variables. The first T constraints are quadratic in * T + 1 variables, and the objective function is quadratic in * T + 1 other variables. * * The control problem models a spring, mass and damper system. * It is of the form * * -------------------------------------------------------------------- * | minimize 1/2 sum x(t)**2 (t = 0 to T) | * | | * | subject to | * | y(t+1) = y(t) - 0.01 y(t)**2 - 0.004 x(t) + 0.2 u(t) | * | | * | x(t+1) = x(t) + 0.2 y(t), | * | | * | y(t) >= -1, -0.2 <= u(t) <= 0.2, | * | | * | (all for t = 0 to T-1) | * | and | * | y(0) = 0, y(T) = 0, x(0) = 10. | * -------------------------------------------------------------------- * * For large enough T (e.g. T >= 90), the optimal objective value * is about 1186.382. * * This model with T = 100 was used as test problem 5.11 in * B. A. Murtagh and M. A. Saunders (1982), A projected Lagrangian * algorithm and its implementation for sparse nonlinear constraints, * Mathematical Programming Study 16, 84--117. * * 14 Nov 1994: First version of spring.f, derived from manne.f. * 26 Sep 1997: Updated for SNOPT 5.3. * ------------------------------------------------------------------ program spring implicit double precision (a-h,o-z) parameter ( maxT = 2000, $ maxm = 2*maxT, $ maxn = 3*maxT + 2, $ maxne = 7*maxT, $ nName = 1 ) character*8 ProbNm , Names(nName) integer T integer ha(maxne) , hs(maxm+maxn) integer ka(maxn+1) double precision a(maxne) , bl(maxm+maxn), bu(maxm+maxn), $ x(maxm+maxn), pi(maxm) , rc(maxm+maxn) * USER workspace (none required) parameter ( lenru = 1) double precision ru(lenru) parameter ( leniu = 1) integer iu(leniu) parameter ( lencu = 1) character*8 cu(lencu) * SNOPT workspace parameter ( lenrw = 1000000) double precision rw(lenrw) parameter ( leniw = 500000) integer iw(leniw) parameter ( lencw = 500) character*8 cw(lencw) logical byname character*20 lfile external SprCon, SprObj * ------------------------------------------------------------------ * Specify some of the SNOPT files. * iSpecs is the Specs file (0 if none). * iPrint is the Print file (0 if none). * iSumm is the Summary file (0 if none). * nout is an output file used here by the main program. iSpecs = 4 iPrint = 15 iSumm = 6 nout = 6 byname = .true. if ( byname ) then * Unix and DOS systems. Open the Specs and print files. lfile = 'spring.spc' open( iSpecs, file=lfile, status='OLD', err=800 ) lfile = 'spring.out' open( iPrint, file=lfile, status='UNKNOWN', err=800 ) end if * ------------------------------------------------------------------ * Set options to their default values. * ------------------------------------------------------------------ call snInit( iPrint, iSumm, $ cw, lencw, iw, leniw, rw, lenrw ) * ------------------------------------------------------------------ * Read a Specs file. This must include "Problem number T" * for some integer T. * ------------------------------------------------------------------ call snSpec( iSpecs, inform, $ cw, lencw, iw, leniw, rw, lenrw ) if (inform .ge. 2) then write(nout, *) 'iSpecs > 0 but no Specs file found' stop end if * ------------------------------------------------------------------ * The following call fetches T, which defines the number of * nonlinear constraints. * It is specified at runtime in the SPECS file. * ------------------------------------------------------------------ call sngeti( 'Nonlinear constraints', T, $ inform, cw, lencw, iw, leniw, rw, lenrw ) if (T .le. 1 .or. T .gt. maxm/2) then write(nout, *) 'Invalid no. of Nonlinear constraints:', T stop end if * Write T into the problem name. write(ProbNm, '(i8)') T if (T .lt. 100) then ProbNm(1:6) = 'Spring' else if (T .lt. 1000) then ProbNm(1:5) = 'Sprng' else if (T .lt. 10000) then ProbNm(1:4) = 'Spri' else ProbNm(1:3) = 'Spr' end if write(nout, *) 'Spring optimal control problem. T =', T * ------------------------------------------------------------------ * Generate an T-period problem. * ------------------------------------------------------------------ call spdata( T, maxm, maxn, maxne, inform, $ m, n, ne, nnCon, nnObj, nnJac, $ a, ha, ka, bl, bu, hs, x, pi ) if (inform .ge. 1) then write(nout, *) 'Not enough storage to generate a problem ', $ 'with Nonlinear constraints =', T stop end if * ------------------------------------------------------------------ * Specify options that were not set in the Specs file. * i1 and i2 may refer to the Print and Summary file respectively. * Setting them to 0 suppresses printing. * ------------------------------------------------------------------ maxS = T i1 = 0 i2 = 0 call snseti( 'Superbasics Limit ', maxS , i1, i2, inform, $ cw, lencw, iw, leniw, rw, lenrw ) * ------------------------------------------------------------------ * Go for it, using a Cold start. * iobj = 0 means there is no linear objective row in a(*). * Objadd = 0.0 means there is no constant to be added to the * objective. * hs need not be set if a basis file is to be input. * Otherwise, each hs(1:n) should be 0, 1, 2, 3, 4, or 5. * The values are used by the Crash procedure * to choose an initial basis B. * If hs(j) = 0 or 1, column j is eligible for B. * If hs(j) = 2, column j is initially superbasic (not in B). * If hs(j) = 3, column j is eligible for B and is given * preference over columns with hs(j) = 0 or 1. * If hs(j) = 4 or 5, column j is initially nonbasic. * ------------------------------------------------------------------ iObj = 0 ObjAdd = 0.0d+0 call snopt ( 'Cold', m, n, ne, nName, $ nnCon, nnObj, nnJac, $ iObj, ObjAdd, ProbNm, $ SprCon, SprObj, $ a, ha, ka, bl, bu, Names, $ hs, x, pi, rc, $ inform, mincw, miniw, minrw, $ nS, ninf, sinf, obj, $ cu, lencu, iu, leniu, ru, lenru, $ cw, lencw, iw, leniw, rw, lenrw ) if (inform .eq. 42 .or. inform .eq. 43 .or. inform .eq. 44) then write(nout, *) ' ' write(nout, *) 'Estimate of required lencw:', mincw write(nout, *) 'Estimate of required leniw:', miniw write(nout, *) 'Estimate of required lenrw:', minrw go to 900 end if write(nout, *) ' ' write(nout, *) 'SNOPT finished.' write(nout, *) 'inform =', inform write(nout, *) 'nInf =', nInf write(nout, *) 'sInf =', sInf write(nout, *) 'Obj =', Obj if (inform .ge. 20) go to 900 stop * ------------------------------------------------------------------ * Error exit. * ------------------------------------------------------------------ 800 write(nout, 4000) 'Error while opening file', lfile stop 900 write(nout, *) ' ' write(nout, *) 'STOPPING because of error condition' stop 4000 format(/ a, 2x, a ) * end of main program for spring problem end