C言語で降圧チョッパ回路における負荷電圧を数値解として算出する

まず、降圧チョッパ回路の負荷電圧波形を取り出せるプログラムをC言語で記述します。
数値解を求める手法は、Runge-Kutta法を使用します。降圧チョッパ回路では、負荷電圧がスイッチングのDuty比に依存します。よって、PI制御では、このDuty比を制御します。初期のDuty比を0.75とし、150[V]に整定した時点で、PI制御を行いDuty比が0.375つまり、負荷電圧が75[V]となる様に各変数を決定しました。比例ゲイン及び、積分ゲインは適性制御、不足制御、過制御の3つの組み合わせのゲインを特定しました。これらの数値解は、Gnuplotで出力します。

数値解算出

DCDC_PI.c

#include<stdio.h>
#include<math.h>    //ライブラリの選択

#define L 7.7e-6
#define C 1
#define R 2.27e-3
#define end 1.5		//マクロ変数の定義

//各ゲインの値

//適正制御
#define kp 10000
#define ki 0.005

////不足制御
//#define kp 0.0009
//#define ki 0.05

////過制御
//#define kp 4
//#define ki 500


double func_IL(double Vo, double E) {

	return (E - Vo) / L;

}	//　インダクタ電流算出用関数

double func_Vo(double IL, double Vo) {

	return (IL - (Vo / R)) / C;

}	//　負荷電圧算出用関数

double func_U( double Vo,double f,double h) {
	
	return (kp * (75 - Vo)) + (ki * (f + ((75 - Vo) * h)));

}	//　PI制御算出用関数


int main() {

	int i, p, error;
	double E, Vo, Io, IL, t,U;
	double k, k1, k2, k3, k4;
	double e, e1, e2, e3, e4;
	double freq, h;
	double f;
	FILE *fp1;											// 変数の定義

	freq = 500;
	h = 1.0 / freq / 100.0;
	Vo = 0.0;
	Io = 0.0;
	IL = 0.0;
	U = 0.0;
	f = 0.0;
	p = 1;												//　各変数の初期化

	if ((error = fopen_s(&fp1, "DCDC_PI.txt", "w")) != 0) {

		printf("\nThis file can't opened \n\n");

		return 0;

	}									// データププロット用ファイル作成(負荷電圧)

	for (t = 0.0, i = 1; t < end; t = t + h, i++) {		// ルンゲクッタ法による演算

		if (t < 0.1) {									//150Vまで挙げる時のDuty比

			if (i > (100 * p - 25) && i < (100 * p)) {

				E = 0;

			}
			else {

				E = 200;

			}

			f = 0;

		}
		else {											//PI制御にかける際のDuty比

			if (i > (((100 * p )-25 ) + (100 * U)) && i < (100 * p) ){

				E = 0;

			}
			else {

				E = 200;

			}

			f += (75 - Vo) * h;							//積分計算

		}
	
			if (i % 100 == 0) {

				p++;

			}

		k1 = func_IL(Vo, E);
		e1 = func_Vo(IL, Vo);							// 1回目の傾き算出

		k2 = func_IL(Vo + (h / 2)*e1, E);
		e2 = func_Vo(IL + (h / 2)*k1, Vo + (h / 2)*e1);	// 2回目の傾き算出

		k3 = func_IL(Vo + (h / 2)*e2, E);
		e3 = func_Vo(IL + (h / 2)*k2, Vo + (h / 2)*e2);	// 3回目の傾き算出

		k4 = func_IL(Vo + (h*e3), E);
		e4 = func_Vo(IL + (h*k3), Vo + (h*e3));			// 4回目の傾き算出

		k = (k1 + (2.0*k2) + (2.0*k3) + k4) / 6;
		e = (e1 + (2.0*e2) + (2.0*e3) + e4) / 6;		// 数値解導出用の傾き算出

		IL = IL + k*h;									// インダクタ電流算出
		Vo = Vo + e*h;									// 負荷電圧算出

		U = func_U(Vo,f,h);								//	Duty比の偏差算出

		fprintf(fp1, "%lf %lf\n", t, Vo);				//ファイルに負荷電圧を書き込む


	}

	return 0;											//プログラムの終了

}

グラフ出力

グラフのデータは上記のプログラムにおいて、テキストファイルとして出力させます。
出力させたプログラムは、Gnuplotに渡し、波形を出力させます。

PLOT.c

se grid
se xrange [0:1.0]
se yrange [0:170]
se xlabel "s[sec]"
se ylabel "Vr[V]"
plot "DCDC_PI.txt" using 1:2 with lines title "PI"