PID controller slowly oscillates out of control
You have some heuristic approach, as copy from textbook or lessons to algorithm. You don't have any wind-up control, hence your integration and differentiation is not affected by changing Kp. 

The Ziegler-Nichols type PID is little bit different

$$u=K_p (\varepsilon + \dfrac{1}{T_i}\int\varepsilon\cdot dt \ + T_d\dfrac{d\varepsilon}{dt}) $$

You can see, that this ZN controller changes the ki,kd also when changing Kp. Your issue is probably related to integration and absence of anti-windup mechanism.

One way is to make an incremental algorithm. You do this by differentiating the whole PID formula and then integrating back like indefinite integral with supplementary integration constant.

$$u=K_p\int (\dfrac{d\varepsilon}{dt} + \dfrac{1}{T_i}\varepsilon\  + T_d\dfrac{d^2\varepsilon}{dt^2}) dt+C $$

This transaltes to incremental form when discretized:

U(k)= dU + U(k-1)

$$u(k)=K_p(\dfrac{d\varepsilon}{dt}+ \dfrac{1}{T_i}\varepsilon\  + T_d\dfrac{d^2\varepsilon}{dt^2}) + u(k-1)$$

The limiting and anti-windup is then implemented very simple

if u(k) &gt; max
 u(k) = max
elseif u(k) &lt; min
 u(k) = min


This form can be further improved by replacing integration with trapz rule and differentiation with 1st order lag filter. However with this type of PID you are not able to limit just P,I,D parts separately.  

A separate elements of P,I,D is actally an Åström controller. It uses double integrator, the second is so called tracking integrator that subtracts the main integrator when the output is saturated. It's more CPU demanding as it uses square root calculation : \$T_{track}=\sqrt{T_i\cdot\ T_d}\$

EDIT:
Controller is described in book: Digital_Self-tuning_Controllers, Bobal et al.

Industrial PID contrller zn3fd source code. The algorithm is incremental : $$u(k)=u(k-1)+\Delta u(k)$$
$$\Delta u(k)= \Delta P+ \Delta I+\Delta D$$

Transfer function:
$$ G(s) = K_p(1+\dfrac{1}{sT_i}+ \dfrac{sT_d}{1+sT_d/\alpha}) $$

Note : \$\alpha\$ is a filtering factor of D-component it is reccomended to be from 4 to 20, default value is 10. Bigger value means less filtering, lower value means more filter, it makes no sense to use low values.

http://www.mathworks.com/matlabcentral/fx_files/8381/1/content/help/html/STCSL.htm

// industrial PID controller, author: Bobal et al.
// parameters: Kp, Ti, Td, alpha, u_in, u_max
// inputs: Setponit, ProcesVar 
// output: u

double ek, ek1, ek2, uk1, uk2, Setponit, ProcesVar;
double Ts, alpha, u_min, u_max;
double Kpu, Tu, Kp, Ti, Td, Tf, cf, ci, cd;
double q0, q1, q2, p1, p2;

//initialzation, calcualte coefficients
ek1=ek2=uk1=uk2=u=0.0;
Tf = Td/alpha;
cf = Tf/Ts;
ci = Ts/Ti;
cd = Td/Ts;
p1 = -4*cf/(1+2*cf);
p2 = (2*cf-1)/(1+2*cf);
q0 = Kp * (1 + 2*(cf+cd) + (ci/2)*(1+2*cf))/(1+2*cf);
q1 = Kp * (ci/2-4*(cf+cd))/(1+2*cf);
q2 = Kp * (cf*(2-ci) + 2*cd + ci/2 - 1)/(1+2*cf);

// ISR PID algo
ek2 = ek1;
ek1 = ek;
ek = (Setpoint - ProcesVar);
uk2 = uk1;
uk1 = u;

u = q0*ek + q1*ek1 + q2*ek2 - p1*uk1 - p2*uk2;
//limit
if   u&gt;u_max u = u_max;
else if u&lt;u_min u = u_min;


NOTE: This controller may not be used with disabled integration.
You can use float instead of double.

EDIT2:

Computer Control: An Overview It has a java source code for the Åström type controller, this won't suffer if the integration is disabled as the incremental type would. The paper describes the functionality in a very compact form, you have to read it. It has a tracking integrator (reset), filtered D-component (same as Bobal's incremental one) and the output is the sum of P,I,D an excellent example of robust controller.



EDIT 3:
If you want to use only the P-control, then you should bias the output with a working point trough a look up table or constant value. With introducing the integral part this is not needed, as the integrator would bias automatically.







simulate this circuit &ndash; Schematic created using CircuitLab