Electronic – Finding coefficients for a PID controller that minimize the integral of the squared error (ISE)

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 infinity and then find the values for the terms k_1, k_2 and 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_2 and k_3 and k_1(?) 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.