subroutine spdata( T, maxm, maxn, maxne, inform, $ m, n, ne, nnCon, nnObj, nnJac, $ a, ha, ka, bl, bu, hs, x, pi ) implicit double precision (a-h,o-z) integer T integer ha(maxne) , hs(maxn+maxm) integer ka(maxn+1) double precision a(maxne) , bl(maxn+maxm), bu(maxn+maxm), $ x(maxn+maxm), pi(maxm) * ------------------------------------------------------------------ * spdata generates data for the "Spring" optimal control problem. * The constraints take the form * c(x) + A*x - s = 0, * where the Jacobian for c(x) + Ax is stored in a(*), and any * terms coming from c(x) are in the TOP LEFT-HAND CORNER of a(*), * with dimensions nnCon x nnJac. * Note that the right-hand side is zero. * s is a set of slack variables whose bounds contain any constants * that might have formed a right-hand side. * * The objective function is * f(x) + d'x * where d would be row iobj of A (but there is no such row in * this example). f(x) involves only the FIRST nnObj variables. * * On entry, * T is the number of time periods. * maxm, maxn, maxne are upper limits on m, n, ne. * * On exit, * inform is 0 if there is enough storage, 1 otherwise. * m is the number of nonlinear and linear constraints. * n is the number of variables. * ne is the number of nonzeros in a(*). * nnCon is the number of nonlinear constraints (they come first). * nnObj is the number of nonlinear objective variables. * nnJac is the number of nonlinear Jacobian variables. * a is the constraint matrix (Jacobian), stored column-wise. * ha is the list of row indices for each nonzero in a(*). * ka is a set of pointers to the beginning of each column of a. * bl is the lower bounds on x and s. * bu is the upper bounds on x and s. * hs(1:n) is a set of initial states for each x (0,1,2,3,4,5). * x (1:n) is a set of initial values for x. * pi(1:m) is a set of initial values for the dual variables pi. * * 14 Nov 1994: First version of spdata. * ------------------------------------------------------------------ parameter ( zero = 0.0d+0, one = 1.0d+0 ) parameter ( bplus = 1.0d+20, dummy = 0.111111d+0 ) * T defines the dimension of the problem. m = T*2 n = T*3 + 2 nb = n + m nnCon = T nnObj = T*2 + 2 ! y(0:T) and x(0:T) nnJac = T + 1 ! y(0:T) ne = T*7 * Check if there is enough storage. inform = 0 if (m .gt. maxm ) inform = 1 if (n .gt. maxn ) inform = 1 if (ne .gt. maxne) inform = 1 if (inform .gt. 0 ) return * ------------------------------------------------------------------ * Generate columns for y(t), t = 0 to T. * The first T rows are nonlinear, and the next T are linear. * The Jacobian is T x (T+1) upper bidiagonal. * We generate the sparsity pattern here. * We put in dummy numerical values for the nonlinear gradients. * The true non-constant values are computed by funcon. * ------------------------------------------------------------------ j = 0 ! counts the variables ne = 0 ! counts the Jacobian and linear constraint entries do 100, k = 0, T j = j + 1 ka(j) = ne + 1 bl(j) = - one bu(j) = bplus x (j) = - one hs(j) = 0 ! Make the y(t) nonbasic. * There are two Jacobian nonzeros per column, * except in the first and last column. if (k .gt. 0) then ! Aij = 1 ne = ne + 1 ha(ne) = k a(ne) = one end if if (k .lt. T) then ! Aij = .02y - 1 (nonlinear) ne = ne + 1 ha(ne) = k + 1 a(ne) = dummy end if * Below the Jacobian the linear constraints are diagonal. if (k .lt. T) then ne = ne + 1 ha(ne) = T + k + 1 a(ne) = - 0.2d+0 end if 100 continue * ------------------------------------------------------------------ * Generate columns for x(t), t = 0 to T. * They form 0.004*I in the first T rows, * and an upper-bidiagonal in the last T rows. * ------------------------------------------------------------------ do 200, k = 0, T j = j + 1 ka(j) = ne + 1 bl(j) = - bplus bu(j) = bplus x (j) = zero hs(j) = 3 ! Make the x(t) basic. * Part of 0.004*I. if (k .lt. T) then ne = ne + 1 ha(ne) = k + 1 a(ne) = 0.004d+0 end if * The bidiagonal parts have two entries * except in the first and last columns. if (k .gt. 0) then ! Aij = 1 ne = ne + 1 ha(ne) = T + k a(ne) = one end if if (k .lt. T) then ! Aij = - 1 ne = ne + 1 ha(ne) = T + k + 1 a(ne) = - one end if 200 continue * ------------------------------------------------------------------ * Generate columns for u(t), t = 0 to T-1. * They form -0.2I in the first T rows. * ------------------------------------------------------------------ do 300, k = 0, T - 1 j = j + 1 ka(j) = ne + 1 bl(j) = - 0.2d+0 bu(j) = 0.2d+0 x (j) = zero hs(j) = 3 ! Make the u(t) basic. ne = ne + 1 ha(ne) = k + 1 a(ne) = - 0.2d+0 300 continue * ka(*) has one extra element. * Some of the variables are fixed. ka(n+1) = ne + 1 bl(1) = zero ! y(0) = 0 bu(1) = zero bl(T+1) = zero ! y(T) = 0 bu(T+1) = zero bl(T+2) = 10.0d+0 ! x(0) = 10 bu(T+2) = 10.0d+0 * ------------------------------------------------------------------ * Set bounds on the slacks. * We don't need to set initial values and states for slacks * (assuming SNOPT does a cold start). * ------------------------------------------------------------------ do 500, j = n + 1, nb bl(j) = zero bu(j) = zero 500 continue * Initialize pi. * SNOPT requires only pi(1:nnCon) to be initialized. * We initialize all of pi just in case. do 600, i = 1, T pi(i) = zero pi(T+i) = zero 600 continue * end of spdata end