Controlling an animatronics bust

Is the arduino millis() function good source for determining dt?

Its the best you will get without using an outside timer.
This method isn't exactly correct, but it is a good approximation. What you are trying to do is figure out what the error will be in the future based on the current slope of the error. So dt should technically be the time step between now and the next step. If your time steps are more or less uniform this isn't a problem. If they are all over the place then you are gonna have problems.

The million dollar dumb question is: what do I do with output now to determine forward or reverse and speed? Should I just have that calculation done right in this routine?

The PID routine takes care of direction. If you have overshot then the output will be in the opposite direction from before because your error will of the opposite sign from before. You set your output to the output variable of your function

Finally, what variables do I need to keep track of outside of this routine for handling the several motors?

You should make your PID a function in which you pass it the set_position, the current_position, the previous_milli, and the previous_error with previous_milli and previous_error being reference variables so that you can set them in the function.

You should really look at the code in the PID library, it will answer a lot of your questions. One thing I'm not so sure if is why they subtract the D term.

/* Compute() **********************************************************************
 *     This, as they say, is where the magic happens.  this function should be called
 *   every time "void loop()" executes.  the function will decide for itself whether a new
 *   pid Output needs to be computed.  returns true when the output is computed,
 *   false when nothing has been done.
 **********************************************************************************/ 
bool PID::Compute()
{
   if(!inAuto) return false;
   unsigned long now = millis();
   unsigned long timeChange = (now - lastTime);
   if(timeChange>=SampleTime)
   {
      /*Compute all the working error variables*/
	  double input = *myInput;
      double error = *mySetpoint - input;
      ITerm+= (ki * error);
      if(ITerm > outMax) ITerm= outMax;
      else if(ITerm < outMin) ITerm= outMin;
      double dInput = (input - lastInput);
 
      /*Compute PID Output*/
      double output = kp * error + ITerm- kd * dInput;
      
	  if(output > outMax) output = outMax;
      else if(output < outMin) output = outMin;
	  *myOutput = output;
	  
      /*Remember some variables for next time*/
      lastInput = input;
      lastTime = now;
	  return true;
   }
   else return false;
}