Hi guys,
i have a problem with NLopt SLSQP in C,
my compiler throws out a exception after 2 iterations of optimization in my
object function, at my computation of gradient, it said:
Exception thrown: write access violation.
grad was nullptr.
the problem is , i used the same code on matlab , it works fine, but the
results of cost function under the initialized value between matlab and c are
not the same, the gradient as well,
the difference is about 1e-5, and the difference between gradient are huge,
could that caused my exception?
i noticed that the difference from xAll(solution of A*x + B*u = x_dot, it is my
state equation) started at 200th row(xAll should be 2500x5), and the difference
is about 1e-12 to 1e-15
did anyone know why?
attached is my c-code , could anyone help me
#pragma warning (disable : 4146)
// nlopt-Bibliothek
#include"nlopt.h"// nlopt Header-Datei
//c-Standard-Bibliothek
#include <stdio.h>
#include <stdlib.h>
#include<math.h>
#include <string.h>
#include<malloc.h>
//selbst geschriebene Bibliothek
#include "light_matrix.h" //Bibliothek für Matrixberechnung
#include"linspace.h"//wie linspace in MATLAB
#include"one_D_Array.h" //generieren dynamische Speicherplatz für 1
dimensionale Matrix
#include"two_D_Array.h"//generieren dynamische Speicherplatz für 2
dimensionale Matrix
#include"KostenFunktion_linear.h"//berechnen den Wert der Kostenfunktion
#include"Berechnungsmodell.h"
#include"struct.h" //Deklaration der Parameter[
double ObjFunc(unsigned n, const double *OptVec, double *grad, void *data)
{
tag*para = (tag*)data;
// Messwerte der zu identifizierenden Parametern
// messdaten of the identified parameter
printf("\ntrue vector is %f %f %f\n", para->trueOptVec[0],
para->trueOptVec[1], para->trueOptVec[2]);
//Kontinuierliche Zeit-Vektor
// time vektor
double *p_contTimeVec = NULL;
p_contTimeVec = one_D_Array(p_contTimeVec, para->NT);
linspace(para->t_start, para->t_end, para->NT, p_contTimeVec);
//berechnen die Lösung des Zustandsmodeles A*X+B*u = x_dot mit
Messwerte der zu identifizierenden Parametern
//solution of the A*x +B*u =x_dot
double **xAll; // 2500 X 5
xAll = Berechnungsmodell( para->trueOptVec, para,\
p_contTimeVec, 0);
//die 2.te und 3.te Spalte der oberen Matrix(xAll), damit rechnet man
den Wert der Kostenfunktion
// the second and third column from xAll, with them to compute the
gradient
Mat phiDotMMeas_init;
Mat *phiDotMMeas;
phiDotMMeas = &phiDotMMeas_init;
MatCreate(phiDotMMeas, para->NT, 1);
Mat phiAMeas_init;
Mat *phiAMeas;
phiAMeas = &phiAMeas_init;
MatCreate(phiAMeas, para->NT, 1);
for (int i = 0; i < para->NT; i++) {
phiDotMMeas->element[i][0] = xAll[i][1];
phiAMeas->element[i][0] = xAll[i][2];
}
//Pointer xAll freigeben
//Point release
for (int i = 0; i < para->NT; i++) {
free(xAll[i]);
}
free(xAll);
xAll = NULL;
//Gradient berechnen
// compute the gradient
/*wie die Syntax in MATLAB: optFun = @(initOptVec)
KostenFunktion_linear(para, initOptVec, model_input_fun, time, \
phiDotMMeas, phiAMeas);*/
// die Werte der Kostenfunktion je nach Methode-FiniteDifferenzes
// values of cost function according to finte difference
double myfunc, myfunc_g1, myfunc_g11, myfunc_g2, myfunc_g22
,myfunc_g3,myfunc_g33;
//pass = 0 , OptVec ist der zu identifzierende Parameter-Vektor
myfunc = KostenFunktion_linear( OptVec, p_contTimeVec, phiDotMMeas, \
phiAMeas, 0,para);
//pass = 1,Vektor x ist {x[0]+h,x[1],x[2]}
myfunc_g1 = KostenFunktion_linear( OptVec, p_contTimeVec, phiDotMMeas, \
phiAMeas, 1, para);
//pass = 11,Vektor x ist {x[0]-h,x[1],x[2]}
myfunc_g11 = KostenFunktion_linear( OptVec, p_contTimeVec, phiDotMMeas,
\
phiAMeas, 11, para);
//pass = 2,Vektor x ist {x[0],x[1]+h,x[2]}
myfunc_g2 = KostenFunktion_linear( OptVec, p_contTimeVec, phiDotMMeas, \
phiAMeas, 2, para);
//pass = 22,Vektor x ist {x[0],x[1]-h,x[2]}
myfunc_g22 = KostenFunktion_linear( OptVec, p_contTimeVec,
phiDotMMeas, \
phiAMeas, 22, para);
//pass = 3,Vektor x ist {x[0],x[1],x[2]+h}
myfunc_g3 = KostenFunktion_linear( OptVec, p_contTimeVec, phiDotMMeas,
\
phiAMeas, 3, para);
//pass = 33,Vektor x ist {x[0],x[1],x[2]-h}
myfunc_g33 = KostenFunktion_linear( OptVec, p_contTimeVec, phiDotMMeas,
\
phiAMeas, 33, para);
//Gradient berechnen je nach Methode-FiniteDifferenzes
//gradient compute
grad[0] = (myfunc_g1 - myfunc_g11)/(2.0*para->h);//after 2 iterations
of optimization : Exception thrown: write access violation:grad was nullptr.
grad[1] = (myfunc_g2 - myfunc_g22)/(2.0*para->h);
grad[2] = (myfunc_g3 - myfunc_g33)/(2.0*para->h);
printf("\n myfunc is %.10f", myfunc);
printf("\n myfunc_g1 is %.10f myfunc_g2 is %.10f myfunc_g3 is %.10f
\n", myfunc_g1, myfunc_g2, myfunc_g3);
printf("\n myfunc_g11 is %.10f myfunc_g22 is %.10f myfunc_g33 is %.10f
\n", myfunc_g11, myfunc_g22, myfunc_g33);
printf("\nmyfunc_g1- myfunc_g11 is %.10f,myfunc_g2- myfunc_g22 is
%.10f,,myfunc_g3- myfunc_g33 is %.10f,\n", myfunc_g1 - myfunc_g11, myfunc_g2 -
myfunc_g22, myfunc_g3 - myfunc_g33);
printf("\nx[0] is %f OptVec[1] is %f OptVec[2] is %f\n", OptVec[0],
OptVec[1], OptVec[2]);
printf("\ngrad[0] is %f grad[1] is %f grad[2] is %f\n", grad[0],
grad[1], grad[2]);
MatDelete(phiDotMMeas);
MatDelete(phiAMeas);
free(p_contTimeVec);
p_contTimeVec = NULL;
return myfunc;
}
int main() {
//Parameter zuweisen
//parameter assignment
tag*para;
tag paraa;
para = ¶a;
para->NT = 2500;
para->NX = 5;
//physikalisches Modell
para->P1 = 2.0;
para->I1 = 32.0;
para->P2 = 4096.0;
para->u = 203.52;
para->R = 764.0;
para->c = 8938100.0 ;
para->J_MR = 1.18/100.0 + 27.3/10000.0;
para->J_AG = 1562.0;
para->maxTorque = 16.9;
para->J_AG_min = 361.0;
para->cG_min = 89381.0;
para->R_min = 7.64;
// dabei werden A und B durch 10 dividiert, damit kann der Programm mit
wenigen Iterationen die Parameter idntifizieren(in MATLAB)
para->val_B[0] = 0.0;
para->val_B[1] = (para->P1 * para->P2 / para->J_MR)/10.0;
para->val_B[2] = 0.0;
para->val_B[3] = 0.0;
para->val_B[4] = 1/10.0;
para->val_A[0] = 0.0;
para->val_A[1] = 1/10.0;
para->val_A[2] = 0.0;
para->val_A[3] = 0.0;
para->val_A[4] = 0.0;
para->val_A[5] = (-para->c / (para->J_MR*pow(para->u, 2)) - para->I1 /
para->J_MR) / 10.0;
para->val_A[6] = (-para->P1 / para->J_MR) / 10.0;
para->val_A[7] = (para->c / (para->J_MR*para->u) - para->P1 * para->P2
/ para->J_MR) / 10.0;
para->val_A[8] = 0.0;
para->val_A[9] = (para->I1*para->P2 / para->J_MR) / 10.0;
para->val_A[10] = 0.0;
para->val_A[11] = 0.0;
para->val_A[12] = 0;
para->val_A[13] = 1/10.0;
para->val_A[14] = 0.0;
para->val_A[15] = (para->c / (para->u*para->J_AG)) / 10.0;
para->val_A[16] = 0.0;
para->val_A[17] = ((-para->c) / para->J_AG) / 10.0;
para->val_A[18] = ((-para->R) / para->J_AG)/10.0;
para->val_A[19] = 0.0;
para->val_A[20] = 0.0;
para->val_A[21] = 0.0;
para->val_A[22] = -1/10.0;
para->val_A[23] = 0.0;
para->val_A[24] = 0.0;
//trueOptVec
para->trueOptVec[0] = 1562.0;
para->trueOptVec[1] = 8938100.0;
para->trueOptVec[2] = 764.0;
//Zeit
para->t_start = 0;
para->dt_init = 0.001;
para->t_end = (para->NT - 1) * para->dt_init;
//Schrittweite der Methode-FiniteDifferenzes
para->h = 1e-07;
double OptVec[3] = {.0};
double initOptVec[3];
//double ww = rand() % 1 + (double)(rand() % 1000) / 1000.0;
//double xx = rand() % 1 + (double)(rand() % 1000) / 1000.0;
//double hh = rand() % 1 + (double)(rand() % 1000) / 1000.0;
//double random_OptVec[3] = { 0 };
//random_OptVec[0] = ww;
//random_OptVec[1] = xx;
//random_OptVec[2] = hh;
////initiale Werte der zu identifizierenden Parameter, dieselbe Werte
in MATLAB, um den Pragramm zu vergleichen
//for (int i = 0; i < 3; i++) {
// initOptVec[i] = para->trueOptVec[i] * (1 + 2 *
(random_OptVec[i] - 0.5));
//}
initOptVec[0] = 0.000443252921871 *pow(10, 6);
initOptVec[1] = 7.539489040483818 *pow(10, 6);
initOptVec[2] = 0.001399243882489 *pow(10, 6);
for (int i = 0; i < 3; i++) {
OptVec[i] = initOptVec[i];
printf("initial OptVec is \n");
printf(" %f", OptVec[i]);
printf("\n");
}
double f_min = 0;//-10000000;
double tol = 1e-11;
double J_AG_min = 361.0;
double cG_min = 89381.0;
double R_min = 7.64;
double lb[] = { J_AG_min, cG_min, R_min };
double ub[] = { 1e+5, 1e+10, 1e+5 };
nlopt_opt opter = nlopt_create(NLOPT_LD_SLSQP, 3);
nlopt_set_min_objective(opter, ObjFunc, para);
nlopt_set_lower_bounds(opter, lb);
nlopt_set_upper_bounds(opter, ub);
nlopt_set_xtol_rel(opter, tol);
nlopt_set_ftol_abs(opter, tol);
nlopt_set_maxeval(opter,50);
nlopt_result result = nlopt_optimize(opter, OptVec, &f_min);
printf("Maximum utility=%f, OptVec=(%f,%f,%f)\n", f_min, OptVec[0],
OptVec[1],OptVec[2]);
nlopt_destroy(opter);
system("pause");
}
//test-code, kleiner Test für NLopt-SLSQP
//double myfunc(unsigned n, const double *x,double *grad,void *data) {
//
// double h = 1e-8;
//
// grad[0] = (log(x[0] + h) - log(x[0])) / h;
// grad[1] = (log(x[1] + h) - log(x[1])) / h;
//
// //grad[0] = 1.0 / x[0];
// //grad[1] = 1.0 / x[1];
// //printf("\ngrad_a[0] is %10f grad_a[1] is %10f\n", grad_a[0],
grad_a[1]);
// printf("\ngrad[0] is %10f grad[1] is %10f\n", grad[0], grad[1]);
// printf("\nx[0] is %10f x[1] is %10f\n",x[0],x[1]);
//
//
// return log(x[0]) + log(x[1]);
//}
//
//double myconstraint(unsigned n, const double *x, double *grad, void*data) {
// double *a = (double *)data;
// grad[0] = a[0];
// grad[1] = a[1];
// return x[0] * a[0] + x[1] * a[1] - 5; //NLopt always expects
constraints to be of the form myconstraint(x) ⤠0
//}
//
//double myinconstraint(unsigned n, const double *x, double *grad, void *data) {
//
// grad[0] = 1;
// grad[1] = -1;
// return x[0] - x[1];
//}
//
// int main(){
//
//
// double f_max = 10000;
// double tol = 1e-6;
// double p[2] = { 1,2 };
// double x[2] = { 1,1 };
// double lb[2] = { 0.001,0.001 };
// double ub[2] = { 10000,10000 };
// // nlopt_opt local_opt ï¼ member from struct, opter
// nlopt_opt opter = nlopt_create(NLOPT_LD_SLSQP, 2); //Optimizer
definieren,2 Diemensionen; struct: nlopt_opt_s , struct variable: nlopt_opt
// nlopt_set_max_objective(opter, myfunc, NULL);// opter steht fuer
algorithmus
//
// nlopt_set_lower_bounds(opter, lb);
// nlopt_set_upper_bounds(opter, ub);
// nlopt_add_equality_constraint(opter, myconstraint, p, tol);// p
parameters
// nlopt_add_inequality_constraint(opter, myinconstraint, NULL, tol); //
Null no raletive parameter
// nlopt_set_xtol_rel(opter, tol);
// nlopt_set_ftol_abs(opter, tol);
// nlopt_set_force_stop(opter, tol);
//
//
// nlopt_result result = nlopt_optimize(opter, x, &f_max);//?
//
// if (result)
// printf("Maximum utility=%f, x=(%f,%f)\n", f_max, x[0], x[1]);
// // free
// nlopt_destroy(opter);
// return 0;
//}
_______________________________________________
NLopt-discuss mailing list
NLopt-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/nlopt-discuss
