Is there any way to speed up many instance of TravelDistance

As flinty said, most of the time is required to download data. There are a couple of changes that will help that problem: 1) get CityData only for states near cityCenter, and 2) use TravelDistanceList instead of TravelDistance.

If we test fewer cities, the process will be faster. We can limit the number of cities (and make fewer calls to CityData) by selecting the states with boundaries close to cityCenter. Here's how to to select the states.

nearbyStates = Select[
   Join[AdministrativeDivisionData[stateCenter,
     "BorderingStates"], {stateCenter}], 
   QuantityMagnitude@
      GeoDistance[cityCenter, #, DistanceFunction -> "Boundary"] < 
     radius*(1 + tolerance) &];

For St. Louis and radius = 50, we call CityData for two states instead of nine.

Next, we need to get distances more efficiently because network calls are the slowest part of the task. TravelDistance needs 1 network call for each distance it finds. For me, finding the distances ran for more than 15 minutes.

Instead, a call to TravelDistanceList returns the all the distances between every pair of locations in a list. We can get many distances with one network call. I found the task was complete in less than 4 minutes.

However, TravelDistanceList doesn't accept long lists of locations (250 seems to work), so the list of cities must be grouped into "edible" chunks. TravelDistanceList returns distances between pairs of locations (saved as distList), but we need only the odd-numbered results. Combine cityList with the odd-numbered distances, and group each city and its distance as table1.

Here's the code to limit nearbyStates and use TravelDistanceList. I've simplified calls to Interpreter[...], and changed the code for cityList and table1.

cityCenter = Interpreter["City"]["St. Louis"];
stateCenter = 
  cityCenter[EntityProperty["City", "AdministrativeDivision"]];
radius = 50;
tolerance = 0.05;
(*remove states if cityCenter is too far from a state's border*)
nearbyStates = Select[
   Join[AdministrativeDivisionData[stateCenter,
     "BorderingStates"], {stateCenter}], 
   QuantityMagnitude@
      GeoDistance[cityCenter, #, DistanceFunction -> "Boundary"] < 
     radius*(1 + tolerance) &];
(*distanceToCenter isn't needed, but it's useful for checking results*)
distanceToCenter[s_] := 
   QuantityMagnitude[TravelDistance[cityCenter, s]]
cityList = Sort[Flatten[
    CityData[{All, ##}] & @@@ 
     EntityValue[nearbyStates, "CanonicalName"]]];
distList = QuantityMagnitude[
   TravelDistanceList /@ Partition[
     Riffle[ConstantArray[cityCenter, Length[cityList]], cityList],
     UpTo[250]]
   ];
table1 = Partition[
   Riffle[cityList, 
    Flatten[#[[Range[1, Length[#], 2]]] & /@ distList]],
   2];
table2 = Join[{cityCenter}, 
   Select[table1, 
    radius*(1 - tolerance) <= #[[2]] <= radius*(1 + tolerance) &]];
graph1 = GeoListPlot[First /@ Rest[table2], PlotMarkers -> Point,
      GeoLabels -> Automatic, GeoRangePadding -> Scaled[0.5],
      ImageSize -> Medium, GeoBackground -> GeoStyling["StreetMap"]];
graph2 = GeoGraphics[GeoMarker[cityCenter]];
Show[graph1, graph2]