Problem with adding unmodeled dynamics to a system
It is a very interesting paper, the steps are very nicely explained. 
In your case the solution should be:

G = G*[eye(4,6); [zeros(2,4), eye(2)+W_tau*delta]];


The reason is that your uncertainty based unmodeled dynamics is the dynamics of the actuator and therefore does not influence any other input. The error you are getting is because you are trying to multiply 6 input system with 2 output transfer function. Or in your case 3 input (I am not really sure why you append 3 to your system equations, it is not in the paper - if you want to stick with it just make sure you add gain 3 to the appropriate input). 

In your case the input vector is:
$$ [w(t),\quad u(t)]' = [d_1(t),\quad d_2(t),\quad \theta_{w1}(t),\quad \theta_{w2}(t),\quad u_1(t),\quad u_2(t)]' $$ 

So your delay transfer function should influence only last two input states \$u_1\$ and \$u_2\$. Therefore you will multiply the system G with:

[eye(4,6); [zeros(2,4), eye(2)+W_tau*delta]]


The same approach you can use later on for the loop-shaping and frequency dependent performance weights \$W_{px}\$ and \$W_{pu}\$. Just remember when you are adding the transfer function G_new to the input of plant P you must multiply P*G_new. And when you are adding the transfer function G_new to the output you need to multiply G_new*P.