It seems that Stata margins command does allow to include a lagged DV as a control. In the following example, which I made by using abdata.dta for this question, marginsplot does not work with model 2 (a lagged dependent variable model), while it works as I expected after running model 1.

My hypotheses are as follows:
H1: Industry demand is positively associated with firm employment (The higher the demand in the present year (t), the higher the employment in the following year (t+1)).
H2: Firm's gross capital strengthens the positive relationship between industry demand and firm employment.

I am attaching do file.

webuseabdata, clear

keep id year n w k ys
order id year

label variable n "firm employment (number of employees) (log(emp)"
label variable w "firm’s wage level (log(wage))"
label variable k "firm’s gross capital (log(cap))"
label variable ys "aggregate output in the firm’s sector (a proxy for demand) (log(indoutpt))"

xtset id year

xtreg F.n n ys k w, fe vce(robust) // testing H1
xtreg F.n n c.ys##c.k w, fe vce(robust) // testing H2

sumy ys k
/* Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
ys | 1031 4.638015 .0939611 4.464758 4.85488
k | 1031 -.4415775 1.514132 -4.431217 3.852441 */

//model 1
xtreg n c.ys##c.k w, fe vce(robust)
margins, at(ys=(4.4(0.2)4.9) k=(-4.4(3)3.9))
marginsplot

//model 2 (lagged dependent variable model)
xtreg F.n n c.ys##c.k w, fe vce(robust)
margins, at(ys=(4.4(0.2)4.9) k=(-4.4(3)3.9))
marginsplot
//error message: "default prediction is a function of possibly stochastic quantities other than e(b)"