What is the straightest line I can make using a linear combination of time series I have 3 processes which generate an output in the form of a time series.

I want to choose a linear combination of the processes that will result in the straightest line possible (I think this property might be called Sinuosity)
A trivial solution would be 0*A + 0*B + 0*C which will always give a straight line,  but this isn't very useful.
How can I calculate the non-trivial solution(s)?
It feels like a Linear Regression problem but I can't see a way to reduce to that.
 A: Proximity of a function $f$ to "straightness" measures the extent to which $f$ is "close" to a linear function of time.  Great flexibility and power to specify straightness can be achieved by extending the linear functions, which are linear combinations of the constant function $1$ and the identity function $t\to t$, to a basis $E$ of the (Hilbert) space of $L^2$ integrable functions of the set of times.
Conventional extensions include the polynomials
$$E = (e_0, e_1, e_2, \ldots, e_k, \ldots) = (1, t, t^2, \ldots, t^k, \ldots)$$
but can be any set of linearly independent functions.  Intuitively, the further out we go into one of these bases, the more we "depart from linearity."
Given $p$ time series 
$$\mathrm{x}_i = ((t_1, x_{i1}), (t_2, x_{i2}), \ldots, (t_j, x_{ij}), \ldots, (t_n, x_{in}))$$
let us compute their projections onto the first $p+1$ elements of this basis using ordinary least squares, giving
$$x_{ij} = b_{i0} + b_{i1}t_j + b_{i2}e_2(t_j) + \cdots + b_{ip}e_p(t_j) + \varepsilon_{ij}.$$
We seek a linear combination of the $\mathrm{x}_i$, with coefficients $\mathbf{\lambda} = (\lambda_1, \lambda_2, \ldots, \lambda_p)$ that is "straightest" in the sense that all the coefficients of $e_j$ for $j\gt 1$ vanish.  That is,
$$\sum_{i=1}^p \lambda_i b_{ij} = 0, \ j = 2, 3, \ldots, p.$$
This is the most we can hope for with $p$ series: if we tried to make one more coefficient vanish, we would have $p$ simultaneous linear equations governing the $\lambda_i$ and usually only $\lambda_i=0$ would be the solution.  By invoking only $p-1$ equations, we are guaranteed to have a system with a nontrivial kernel.
We are left to choose a basis element of that kernel.  Assuming it is just one-dimensional (which will generically be the case), we may impose one more condition.  A convenient one is to make the resulting linear combination look like an arithmetic mean of the time series.  I do that by standardizing it so that its variance is $1/p$ times the collective variance of all the time series data.  This can easily be done in two steps: first, require that the coefficients sum to unity:
$$\lambda_1 + \lambda_2 + \ldots + \lambda_p = 1.$$
Then, perform the standardization.  All other solutions will be linear functions of this one plus a linear function of time.
Let's turn to a worked example.  This proposal is implemented in the following R code, which generates an array of time series (which may have irregular spacing and multiple observations per time), finds the coefficients $b_{ik}$ using least squares fits, adjoins the sum-to-unity vector $(1,1,\ldots, 1)$ to this matrix, and solves for $\lambda$.  That directly produces the linear combination of the original series, $\sum_{i}\lambda_i \mathrm{x}_i$, which is then standardized.  The original series are plotted (in color, using dashed lines) and the "straightest" linear combination is overplotted (black solid line) for comparison.

#
# Specify the problem.
#
n <- 96            # Number of time steps
k <- 1             # Observations per time step
p <- 5             # Number of series
shape <- 3         # Higher (positive) values make times more evenly spaced
set.seed(17)       # Makes the results reproducible
#
# Create time series, one per column of `y`.  The times themselves are in `times`.
#
q <- p
times <- rep(c(cumsum(rgamma(n, shape, shape)), NA), k)
n.k <- length(times)
beta <- round(matrix(rnorm(q*p), q), 1)
y <- matrix(rnorm(q*n.k), ncol=p) %*% beta + 
  100 * outer(times, 1:p, function(i,j) sin(2*j*i/n)) + rnorm(p*n.k, sd=sqrt(n.k))
#
# Construct fits using a basis of orthogonal polynomials.
#
x <- times; x[is.na(times)] <- 0 # (`poly` chokes on NA values, so zero them out)
basis <- poly(x, degree=p-1)     # Includes a linear term
fits <- apply(y, 2, function(z) coef(lm(z ~ basis)))
fits["(Intercept)", ] <- 1       # Make coefficients sum to unity
lambda <- solve(fits, c(1, rep(0, p-1)))
y.hat <- y %*% lambda
#
# Standardize the combination to look like an average of the time series
#
y.hat <- (y.hat - mean(y.hat, na.rm=TRUE)) / sd(y.hat, na.rm=TRUE) * 
  sd(y, na.rm=TRUE) / sqrt(p) + mean(y, na.rm=TRUE)
#
# Display the results.
#
colors <- as.list(rainbow(p))
times.range <- range(times, na.rm=TRUE)
plot(times.range, range(c(y, y.hat), na.rm=TRUE), type="n",
     xlab="Time", ylab="Value", main="Data and Fit")
invisible(mapply(function(z, c) lines(times, z, col=c, lwd=2, lty=3), 
                 as.data.frame(y), colors))
lines(times, y.hat, lwd=2, col="Black", lty=1)

A: In a fairly practical sense, one nice and simple way to get "a linear combination close to a straight line" in one particular sense -- at least from data sampled at regular intervals -- would be to regress the time index (1,2,...) on the series A, B, C. That is, fit a least squares regression y ~ A + B + C where y is simply a linearly increasing set of values.
It will be as close as it can get (in the least squares sense) to the chosen line, but it won't necessarily be especially smooth.
You can then scale the coefficients of A,B, and C by a common arbitrary constant, or shift the intercept (e.g. to zero, if you like). 
(I make no claims that doing what you're attempting is meaningful, but if that's really what you want to do, this straight OLS regression approach would be the first thing I'd try - it's simple to do, and easy to explain to people. If you want to minimize something other than the plain sum of squares of deviations from a straight line, then of course there are ways to do that.)

Example: I start with whuber's matrix y from his code.
> ys=as.data.frame(cbind(y,1:97))
> colnames(ys)<-c(LETTERS[1:5],"t")
> tfit=lm(t~.,ys)
> plot(fitted(tfit),type="l")



