If this had to be the approach to measure the impulse response, here are my further suggestions specific to the qeustions:
How to determine the appropriate excitation signal duration?
The equivalent noise bandwidth of each FFT bin if not windowed is $1/T$ where $T$ is the duration of the signal in time prior to taking the FFT (note zero padding will not change this, but windowing will). So this will be the equivalent to the resolution bandwidth on network and spectrum analyzers, and will effect the resolution of the resulting transfer function in frequency. Resolution is not equivalent to number of samples, since we can zero pad in time which will result in more samples in frequency (over the same frequency duration) but this will just be a smoothing or interpolation of the same frequency domain result and will not help us further distinguish a closely spaced variation in the transfer function that may actually exist. To be able to measure such a tight variation in frequency, we need a longer time sample. Ultimately this is the expected length of the impulse response.
Do you recommend excitation signal windowing?
This depends on the excitation signal and the dynamic range of the resulting frequency transfer function. It would generally be a good idea but significantly complicates the processing; for example if you did end up using a swept sine as the excitation signal, the window will be significantly attenuating all frequencies represented at the start and end of the sweep, thus reducing the SNR and therefore quality of the results in the frequency domain transfer function for those frequencies. Windowing decreases the frequency resolution of the main lobe of each DFT bin's equivalent filter at the benefit of significantly decreasing the sidelobes and thus improving overall dynamic range. Please see this post further detailing windowing considerations as this is what would all apply here:
How to calculate resolution of DFT with Hamming/Hann window?
Why would one use a Hann or Bartlett window?
How to determine appropriate output signal measurement duration?
Before applying fft, how much zero padding do I need to add to
excitation signal and output signal?
Both these questions affect the zero padding since we need the FFT's to be of the same duration prior to computing the frequency domain transfer function. The output signal will be the duration of the input signal plus the length of the impulse response. We already established that the input signal should be the duration of the input response, so in this case it would be zero padded out an additional amount. There is no harm in having a longer duration sequence, these are the minimums so if I were to do this it would likely be an iterative approach where I would evaluate both the sampling rate and signal duration from the frequency domain transfer function and time domain impulse response respectively: In the frequency domain, if the transfer function has not sufficiently attenuated by the time we reach half the sampling rate, aliasing will occur and the sampling rate should therefore be increased. In the time domain if the impulse response has not sufficiently attenuated by the time we reach the end of the duration of the resulting answer, time domain aliasing will also occur and the overall time duration of the signal should be increased. To summarize with consideration to signal duration and zero padding: The signal duration must be longer in time than the impulse response duration. This signal must also be zero padded out further by the duration of the impulse response duration to allow for the output to fully settle. The output duration will match the input duration and no further zero padding will be required. The sampling rate must also be confirmed to be sufficiently high enough by reviewing the resulting frequency transfer function.
How to avoid zero division when excitation signal max frequency <
f_s/2 ?
Since you are creating the excitation signal you will know precisely if any frequency is 0, so this would be trivial to avoid. What I would be tempted to do personally is decide what dynamic range my solution needs to cover for the resulting frequency domain transfer function and then replace any frequency domain zeros with a very small number that I know is much smaller than that. However having an excitation signal that had lower energy in any given bin would be a poor choice for determining the frequency transfer function, which is required to be accurate to determine the impulse response. Granted in the design of this you would want to have an excitation signal with a bandwidth that sufficiently exceeds the bandwidth of the transfer function, in which case we could have its energy taper off at the higher frequencies that we know are well beyond the transfer function's bandwidth--and then it that case we would simply ignore (zero) the output beyond that frequency where the result is meaningless (and we wouldn't want it to be applied to the IFFT).
OTHER APPROACHES
There are a lot of drawbacks and challenges to determining an impulse response using this approach similar to the pitfalls in designing filters using the frequency sampling approach detailed here:
Difference between frequency sampling and windowing method
Alternate methods that would be more straightforward would be to slowly sweep a frequency tone in the time domain with a rate less than the resolution bandwidth of the measurement desired while comparing the amplitude and phase of the output and input directly in time (as done in a network analyzer). Alternatively but similarly to step through each frequency and compare the static results for each one. Another approach would be to use a PRN sequence as the excitation tone or any other pseudo-random white noise source (random number generator) and use the Wiener-Hopf equations to determine the frequency domain transfer function which results in the least-mean-square solution for the transfer function. This is further detailed at this post with practical MATLAB code specific to doing this:
How determine the delay in my signal practically
