The problem

The problem is to find the coefficients for the P, I and D terms of a PID controller used to regulate the object that minimize the the integral of the squared error (ISE):

$$Q = \int_0^\infty \epsilon^2(t) \ dt  $$

The object being controlled is described by the transfer function:
$$ G(s) = \frac{1}{s^3 + 6s^2 + 5s} $$

My attempt

My idea was to find to find the expression for error in the s-domain, then perform an inverse Laplace transform to get the expression for the error in the time domain, then integrate the error from $0$ to $\infty$ and then find the values for the terms $K_1, K_2, K_3$ that will minimize the function. I apologize for the mistakes if I've made any as my knowledge of control theory is very basic.

What I managed to do:

The transfer function of a PID controller is:

$$ K(s) = k_1 + \frac{k_2}{s} + k_3 s $$

If $H(s)$ is the overall transfer function, $W(s)$ is the input, let's say it's the unit step function, the expression for error in the s-domain is, if I'm correct:

$$E(s) = W(s) - H(s)W(s)$$
$$E(s) = W(s)[1 - H(s)]$$
$$E(S) = \frac{1}{s}[1 - H(s)]$$

The overall transfer function will be given by:

$$ H(s) = \frac{K(s)G(s)}{1 + K(s)G(s)} $$

After substitution and simplification:

$$ H(s) = \frac{k_3 s^2 + k_1 s + k_2}{k_3 s^2 + k_1 s + k_2 + s^4 + 6 s^3 + 5 s^2} $$

Then, substituting for the expression for error:

$$E(S) = \frac{1}{s}[1 - \frac{k_3 s^2 + k_1 s + k_2}{k_3 s^2 + k_1 s + k_2 + s^4 + 6 s^3 + 5 s^2}]$$

I used Matlab ilaplace command to calculate the inverse Laplace transform of $E(s)$ and got:

$$ \epsilon(t) = \frac{\delta(t)}{t} - \frac{(exp[(- k_2 - t^2 k_3 - 5t^2 - 6t^3 - t^4)(t^2 + 6t + 5)]}{t} $$

I don't know if I've arrived at the correct result, what seems weird to me is that the resulting expression for $\epsilon(t)$ doesn't contain any $k_1$ term. But assuming this is correct, I would finish by finding the values of parameters $k_1$, $k_2$ and $k_3$ that minimize this integral: $Q = \int_0^\infty \epsilon^2(t) \ dt  $ 

Question

Now, I would be grateful if you could verify whether my approach and solution is correct and describe how you would find the values of the parameters that minimize this integral using Matlab i.e. what commands and how you would use them. I'm wondering if lack of closed-form solution to the integral would make finding the parameters impossible.

Or maybe, there is a simpler method of solving this problem, perhaps by avoiding computing the integral.

Those are counterweights. They work exactly the same as those lead counterbalance weights on the wheels of your automobile. 

If they left those out, then those connecting rods and bearings would create an out-of-balance condition, and the wheels would vibrate at higher speeds. That could very well damage the wheels. But as a couple of others have nicely pointed out, the violent shaking could derail the train.
