And similarly, why we have $$\langle x|P|x' \rangle =-i\hbar \frac{d}{dx} \delta(x-x')~?$$
Both equations come from the Wikipedia page Momentum operator.
My question is, if we insert the identity operator,wouldn't it be $$P=\int P |x\rangle\langle x|dx ~?$$ And in the previous equation, seems we just let $$\langle x|P|x' \rangle =P\langle x|x' \rangle=-i\hbar \frac{d}{dx} \delta(x-x').$$ How can we change the position of operator in the expression arbitrarily?
The momentum operator is usually defined by its matrix elemenets, $\left\langle x\right|P\left|\psi\right\rangle :=-i\hbar\frac{d\psi\left(x\right)}{dx}$ where $\left|x\right\rangle$ is a position basis ket, $\left|\psi\right\rangle$ is any state, and $\psi\left(x\right):=\left\langle x|\psi\right\rangle$.
The reason is, one can show the operator $P$ defined in this way is Hermitian, and it satisfies the classical relation $\left\langle \psi\left(t\right)\left|P\right|\psi\left(t\right)\right\rangle =m\frac{d}{dt}\left\langle \psi\left(t\right)\left|X\right|\psi\left(t\right)\right\rangle$ (the proof is via Ehrenfest's theorem using the Hamiltonian $H=\frac{P^{2}}{2m}+V\left(X\right))$.
The position kets $\left|x\right\rangle$ are orthogonal (as they are eigenstates of a Hermitian operator, $X$, all with distinct eigenvalues), and they are normalized such that $\left\langle x|x'\right\rangle =\delta\left(x-x'\right)$ (this is the continuous version of orthonormality).
Then if we use the definition of $P$ for the case of $\left|\psi\right\rangle =\left|x'\right\rangle$, we get $\boxed{\left\langle x\right|P\left|x'\right\rangle =-i\hbar\frac{d}{dx}\delta\left(x-x'\right)}$. It is perfectly okay to differentiate the delta function, in fact you can differentiate it arbitrarily many times (see e.g. https://en.wikipedia.org/wiki/Dirac_delta_function#Derivatives).
Now since $P$ is Hermitian, it has an orthogonal basis of eigenstates $\left|p\right\rangle$ (just like the position basis). They satisfy the eigenvalue equation $P\left|p\right\rangle =p\left|p\right\rangle$ and are normalized such that $\left\langle p|p'\right\rangle =\delta\left(p-p'\right)$. We make use of the resolution of the identity, which is a property of any orthonormal basis, and for the momentum basis it is $\mathbb{I}=\int_{-\infty}^{\infty}\left|p\right\rangle \left\langle p\right|dp$. Therefore $P=P\int_{-\infty}^{\infty}\left|p\right\rangle \left\langle p\right|dp$. The linear operator $P$ is acting on a sum (integral), so let's move it inside the sum $P=\int_{-\infty}^{\infty}P\left|p\right\rangle \left\langle p\right|dp$ and now using the eigenstate property, $\boxed{P=\int_{-\infty}^{\infty}p\left|p\right\rangle \left\langle p\right|dp}$. Any Hermitian operator can be written in this way in terms of its eigenvalues and eigenstates, e.g. $X=\int_{-\infty}^{\infty}x\left|x\right\rangle \left\langle x\right|dx$. By using the resolution of the identity again but in the position basis, the relation $\boxed{P=\int_{-\infty}^{\infty}P\left|x\right\rangle \left\langle x\right|dx}$ holds as well.
Now for the last relation, $\left(-i\hbar\int_{-\infty}^{\infty}\left|x\right\rangle \frac{d}{dx}\left\langle x\right|dx\right)\left|\psi\right\rangle =-i\hbar\int_{-\infty}^{\infty}\left|x\right\rangle \frac{d}{dx}\psi\left(x\right)dx=\int_{-\infty}^{\infty}\left|x\right\rangle \left(-i\hbar\frac{d}{dx}\psi\left(x\right)\right)dx=\int_{-\infty}^{\infty}\left|x\right\rangle \left\langle x\left|P\right|\psi\right\rangle dx=P\left|\psi\right\rangle$
This is true for any $\left|\psi\right\rangle$, so $\boxed{P=-i\hbar\int_{-\infty}^{\infty}\left|x\right\rangle \frac{d}{dx}\left\langle x\right|dx}$. We're not really differentiating a bra, but rather this derivative bit is defined by the way it acts on any state: $\left(\frac{d}{dx}\left\langle x\right|\right)\left|\psi\right\rangle =\frac{d}{dx}\left\langle x|\psi\right\rangle =\frac{d\psi\left(x\right)}{dx}$, so it's easier to make sense of this relation by having it act on a state.
As you assert, we cannot change the order of kets, bras and operators: $\left\langle x\left|P\right|x'\right\rangle \neq P\left\langle x|x'\right\rangle$. On the left hand side, we have an inner product between $\left\langle x\right|$ and $P\left|x'\right\rangle$ , resulting in a scalar. On the right hand side, we have $P$ multiplied by a scalar, resulting in an operator. We can however move scalars around freely, e.g. $cX=Xc$ and $c\left|\psi\right\rangle =\left|\psi\right\rangle c$.
The other answers and comments have made it clear that you can define momentum operator to satisfy $$\tag1\left<x|\hat p|\psi\right>=-i\hslash\frac{\partial\ }{\partial x}\left<x|\psi\right>$$
My contribution will only be that it is actually down to choice. You get to choose between one of many possible definitions. However, they are all equivalent in the end.
One popular alternative is to choose $$\tag2\left<x|p\right>=\alpha e^{ipx/\hslash}$$ for some $\alpha$, be it $1,\frac1{\sqrt L},\frac1{\sqrt{2\pi\hslash}}$ or whatever else is the convention chosen. In my case, I always choose $\alpha=1$ and for the momentum differential element to appear as the combination $\frac{\mathrm dp}{2\pi\hslash}$, so, for example, the identity operator is $\int\left|p\right>\frac{\mathrm dp}{2\pi\hslash}\left<p\right|=\hat{\mathbb I}=\int\left|x\right>\mathrm dx\left<x\right|$
When you pick this wavefunction definition as the postulate, then, by sandwiching the momentum operator $\hat p=\int\left|p\right>\dfrac{p\,\mathrm dp}{2\pi\hslash}\left<p\right|$ between two position representations of the identity operator as shown above, you can derive Equation (1) as required.
In the famous textbook by Feynman and Hibbs, they started only with path integrals in position space. By considering momentum as mass multiplied by velocities, i.e. by taking time data and detector position data, and taking the ratio, they define the measurement of velocities, and thereby the measurement of momentum. Doing that, they obtain Equation (2) and proceed from there. So, in a sense, you do not need to assume either of Equation (1) or (2), you could just work with path integrals and positions only.
From the PoV of Noether's theorem, we want the momentum operator to be the position translation of the wavefunction. i.e. the translation operator must be an exponential of the momentum operator, and the momentum operator must be Hermitian. This then pins it down to within a position function (chosen as zero) and a choice of sign. Note that by this I mean that Taylor's expansion theorem already fixes the translation operator, and Noether's theorem fixes that it must be related to momentum operator.
From the PoV of the rest of quantum mechanics, you kind of want to respect the Stone-von-Neumann theorem. After all, it is the one thing that pretty much fixes the Hilbert space and is a wonderful godsend that is so sorely lacking in QFT.
Between this and Canonical Quantisation, you have $[\hat x{\,}^j,\hat p_k]=i\hslash\,\delta^j_k$ and $\hat x=x$ then pins down the momentum operator up to a function of position. This function of position is usually chosen as zero, but in the radial momentum case, Hermitian-ness requires that that function of position be something else. Once you know it is non-zero, it is easy to go look it up and find out what it has to become. Dirac made that quite clear.
As you can see, there are many inroads to this issue and they converge onto the same mathematical expressions. For some of the above, e.g. Noether's theorem and Stone-von-Neumann, check out the comment linked by Qmechanic.
