How do you turn the output of a nnet neural network model into an equation? 
Assuming the output of the above nnet feedforward model (nnetModel) is such that the following summary is produced:
> summary(nnetModel)
a 2-2-1 network with 9 weights
options were - linear output units decay=0.01
  b->h1  i1->h1  i2->h1
    6      4       5 
  b->h2  i1->h2  i2->h2
    9      7       8 
  b->y   h1->y   h2->y
    3      1       2 

Is there a way of turning this information into an equation which also incorporates the weight decay (0.01)? I would like to make use of this equation as an objective function for nonlinear optimisation.  
I've had a go at converting the model output into an equation function (FFNN) which makes use of the sigmoid activation function (extracted from the C-sources; filennet.c, lines 157-165).
sigmoid <- function(sum){
    if (sum < -15.0) {
        return (0.0) 
    } else if (sum > 15.0) { 
        return (1.0)
    }else 
        return (1.0 / (1.0 + exp(-sum)))
}


FFNN <- function(i1, i2){
    Sumh1 = (i1*4) + (i2*5) + 6
    Sumh2 = (i1*7) + (i2*8) + 9

    HU1 <- sigmoid(Sumh1)
    HU2 <- sigmoid(Sumh2) 

    Sumy = (HU1*1) + (HU2*2) + 3

    Y <-  sigmoid(Sumy)
    print(Y)     
}

I would like to extend this out for my actual nnet model output which has 20+ independent variables.
Does anyone know a better way of writing an equation function which also factors in the weight decay (0.01)?
 A: It would be better if you use vectorization in your code. In that way, FNNN code reduces to simple matrix multiplications and non-liniearity. Here is a simple example in your case.
$ \theta^{(1)}*x = z^{(1)} = \begin{bmatrix}
       \ b_{1}^{(1)} & \ W_{11}^{(1)} & \ W_{12}^{(1)} \\[0.3em]
       \ b_{2}^{(1)} &\ W_{21}^{(1)} & \ W_{22}^{(1)}  
     \end{bmatrix} * \begin{bmatrix}
       \ 1  \\[0.3em] \ x_1  \\[0.3em]
       \ x_2   
     \end{bmatrix} =  \begin{bmatrix} \ z_1^{(1)}  \\[0.3em] \ z_2^{(1)}  \end{bmatrix}$
$ a^{(2)}= sigm(z^{(1)})  $
$ \theta^{(2)}*a^{(2)} = z^{(2)} = \begin{bmatrix}
       \ b_{1}^{(2)} & \ W_{11}^{(2)} & \ W_{12}^{(2)}  
     \end{bmatrix} * \begin{bmatrix}
       \ 1  \\[0.3em]
       \ a^{(2)}_1  \\[0.3em]
       \ a^{(2)}_2   
     \end{bmatrix} = \begin{bmatrix}
       \ z_1^{(2)} 
     \end{bmatrix} $
$ a^{(3)}= sigm(z^{(2)})  $
The feed-forward equation should not include weight decay which is a part of cost function not FF structure. But I don't understand why you want to use the FF equation in another optimization procedure. Isn't it already learned by using gradient decent?
A: Concur with @GeoMatt22 This is exactly what ADs such a Tensorflow, pytorch autograd etc handle for you.  Just write the loss function in tf/autograd, and it will handle backprop, training and so on for you.
Weight decay is as simple as adding a term like weight * weight to your loss function.
