Formula 2) is the correct one in general except it's the derivative of the sectional curvature i.e of $\frac{k_t(X,Y)}{|X\wedge Y|^2_t}$ (and not just of $k_t(X,Y)$) for an orthonormal frame $X,Y$ with respect to the original metric. This accounts for the last term in formula 2.
For $k'(0)$ itself the correct formula is
$$
k'(0)=\nabla_X\nabla_Y h(X,Y)-\tfrac12\nabla_X\nabla_X h(Y,Y)-\tfrac12\nabla_Y\nabla_Y h(X,X)+h(R(X,Y)Y,X)
$$
As Deane said, this is a straightforward calculation. But since the OP seems to be struggling with it I'll supply some details. Let $\nabla^t$ be the Levi-Civita connection for $g_t$. Then $\nabla^t_XY=\nabla_XY+tS(X,Y)+O(t^2)$ where $S$ is a (2,1)-tensor and $\nabla=\nabla^0$.
Let's first compute $k'(0)$ in terms of $S$. As usual let's work in normal coordinates around $p$ with $X,Y$ coordinate fields.
We have $$\langle R^t(X,Y)Y,X\rangle_t=\langle R^t(X,Y)Y,X\rangle_0+t\cdot h(R(X,Y)Y,X)+O(t^2)=$$
$$=\langle \nabla^t_X\nabla^t_YY,X\rangle_0-\langle \nabla^t_Y\nabla^t_XY,X\rangle_0+t\cdot h(R(X,Y)Y,X)+O(t^2)$$
Next we expand the first term
$$\langle \nabla^t_X\nabla^t_YY,X\rangle_0=\langle \nabla^t_X(\nabla_YY+tS(Y,Y)),X\rangle_0+O(t^2)=$$
$$\langle \nabla_X\nabla_YY,X\rangle_0+t\langle S(X,\nabla_YY)+\nabla_XS(Y,Y), X\rangle_0+O(t^2)=$$
$$ \langle \nabla_X\nabla_YY,X\rangle_0+t\langle \nabla_XS(Y,Y), X\rangle_0+O(t^2)$$
where in the last equality we used that $\nabla_YY(p)=0$.
After a similar computation for $\langle \nabla^t_Y\nabla^t_XY,X\rangle_0$ we get that
$$k'(0)=\langle\nabla_XS(Y,Y)-\nabla_YS(X,Y), X\rangle_0+h(R(X,Y)Y,X)$$
$$=\nabla_X S(Y,Y,X)-\nabla_YS(X,Y,X)+h(R(X,Y)Y,X)$$
where we lowered the index and turned $S$ into a $(3,0)$-tensor $S(X,Y,Z)=\langle S(X,Y),Z\rangle_0$
Lastly, recall that
$\langle \nabla_XY,Z\rangle=\frac{1}{2}[X\langle Y,Z\rangle+Y\langle X,Z\rangle -Z\langle X, Y\rangle]$ for coordinate fileds. This easily gives
$S(Y,Y,X)=\frac{1}{2}[Yh(X,Y)+Yh(X,Y)-Xh(Y,Y)]=\nabla_Yh(X,Y)-\frac{1}{2}\nabla_Xh(Y,Y)$ and
$S(X,Y,X)=\frac{1}{2}\nabla_Yh(X,X)$ written invariantly as tensors.
Altogether this gives
$$
k'(0)=\nabla_X\nabla_Y h(X,Y)-\tfrac12\nabla_X\nabla_X h(Y,Y)-\tfrac12\nabla_Y\nabla_Y h(X,X)+h(R(X,Y)Y,X)
$$ as promised.
Formulas 1) and 3) are not true in general but might be true in the specific circumstances where they are applied in the papers in question.
The comment of Robert Bryant essentially answered the question; all the relevant information can be found in his other answer: https://mathoverflow.net/a/100372/39348. To avoid leaving the question open, I reproduce it here.
In all that follows, $T$ and $T^*$ are shorthand for $T_x M$ and $T_x M^*$.
First, a small correction: what I really want to define is a quadratic form on the space of decomposable tensors, that we shall denote by
$$ \Pi_2 := \{X \wedge Y \in \Lambda^2T \;\;|\;\; X, Y \in T\}, $$
rather than on the Grassmannian $\operatorname{Gr}(2, T)$ (which is the projective space corresponding to $\Pi_2$.)
Such an object does makes sense; the appropriate language here is that of algebraic geometry. Note that a tensor $\alpha$ belongs to $\Pi_2$ iff we have $\alpha \wedge \alpha = 0$; this gives a family of quadratic equations that define $\Pi_2$. But there is a classical definition of the ring of polynomials defined on an algebraic variety, as the quotient of the general polynomial ring by the ideal of equations satisfied by the variety. Thus if we denote by $Q(E)$ the space of all quadratic forms on a vector space $E$, then it is sensible to define
$$ Q(\Pi_2) := Q(\Lambda^2 T) / I, $$
where $I$ is the set of elements of $Q(\Lambda^2 T)$ that are identically zero on $\Pi_2$. By the previous remark, $I$ is exactly the image of the map $\Phi: \Lambda^4 T^* \to S^2(\Lambda^2 T^*)$ (the polarization identity allows us to identify this codomain with $Q(\Lambda^2 T)$), defined by the identity
$$ \Phi(\phi)(A, B) := \phi(A \wedge B) $$
for every $\phi \in \Lambda^4 T^*$ (seen as a linear form on $\Lambda^4 T$) and every $A, B \in \Lambda^2 T$. Thus we have
$$ Q(\Pi_2) = Q(\Lambda^2 T) / \operatorname{Im} \Phi. $$
But the map $\Phi: \Lambda^4 T^* \to S^2(\Lambda^2 T^*)$ has a natural left inverse $\Psi: S^2(\Lambda^2 T^*) \to \Lambda^4 T^*$, defined by the identity
$$ \Psi(q)(X \wedge Y \wedge Z \wedge W) := \frac{1}{3} \left( q(X \wedge Y, Z \wedge W) + q(X \wedge Z, W \wedge Y) + q(X \wedge W, Y \wedge Z) \right) $$
for any symmetric bilinear form $q \in S^2(\Lambda^2 T^*)$ and any vectors $X,Y,Z,W \in T$. (We could also define it by $\Psi(\alpha \otimes \beta) = \alpha \wedge \beta$ for any $\alpha, \beta \in \Lambda^2 T^*$).
It follows that $Q(\Pi_2)$ is canonically isomorphic to $\operatorname{Ker} \Psi$. But the condition $\Psi(q) = 0$ is exactly the Bianchi identity.
Best Answer
For the first question: positive curvature operator on a compact manifold implies that the manifold is diffeomorphic to a space form, i.e., a manifold of sectional curvature one. This is due to C. Boehm and B. Wilking Manifolds with positive curvature operators are space forms Annals of Mathematics, 167 (2008), 1079–1097.
Thus most known positively curved manifolds do not admit metrics with positive curvature operator. The simplest example is $CP^n$, $n>1$.
I do not know know much about how sectional curvature pinching can imply positive curvature operator but for example, you could trace references from J.-P. Bourguignon, H. Karcher, Curvature operators: pinching estimates and geometric examples in Ann. Sci. École Norm. Sup. (4) 11 (1978), no. 1, 71–92 where some estimates on the eigenvalues of the curvature operator in terms of sectional curvature pinching can be found.