Adding current files

Adding more files to repository

Easy/Main.m

@@ -0,0 +1,62 @@
+function [X,F,INFO] = Main()
+% clear all
+% Simple snopt application
+global problem
+x0 = 100;
+
+% Variables
+problem.N = 10; % Nodes
+problem.ndof = 5;
+problem.nstates = 2*problem.ndof;
+problem.ncontrols = 3;
+problem.nvarpernode = problem.nstates + problem.ncontrols;		% number of unknowns per time node
+problem.nvar = problem.nvarpernode * problem.N;		%Number of variables (controls and states per time node)
+problem.ncon = problem.nstates * problem.N;
+
+X0 = x0*randn(problem.nvar,1);
+L = -inf*ones(problem.nvar,1);
+U = inf*ones(problem.nvar,1);
+
+FL = [-inf;zeros(problem.ncon,1)];
+FU = [inf;zeros(problem.ncon,1)];
+userfun = 'easyobj';
+
+problem.easyobj = @easyobj;
+% Output informative files
+snprint   ('probName.out');
+snsummary ('prName.sum');
+[X,F,INFO] = snopt(X0,L,U,FL,FU,userfun);
+snprint    off;
+snsummary off;
+
+end
+
+function [F,g]= easyobj(x)
+
+global problem
+N = problem.N;
+
+h = 0.01;
+
+g = zeros(1+problem.ncon,problem.nvar);
+
+F(1) = sum(x(problem.ncon+1:end).^2);
+g(1,:) = 2*[zeros(1,problem.ncon), transpose(x(problem.nstates*N+1:end))];
+
+for i = 1:N
+    if i < N
+        x1 = [x(problem.ndof*(i-1)+(1:problem.ndof));x(problem.ndof*(N+(i-1))+(1:problem.ndof))];
+        x2 = [x(problem.ndof*i+(1:problem.ndof));x(problem.ndof*(N+i)+(1:problem.ndof))];
+    else
+        x1 = [x(problem.ndof*(i-1)+(1:problem.ndof));x(problem.ndof*(N+(i-1))+(1:problem.ndof))];
+        x2 = [x(1:problem.ndof);x(problem.ndof*N+(1:problem.ndof))];
+    end
+    [f, dfdx, dfdxdot, dfdu] = easydyn1(x2,(x1-x2)/h); %(x1+x2)/2
+
+    F = [F;f];
+
+    g(problem.nstates*(i-1)+1+(1:problem.nstates),(i-1)+(1:problem.nstates)) = -dfdxdot/h;
+    g(problem.nstates*(i-1)+1+(1:problem.nstates),i+(1:problem.nstates)) = dfdx + dfdxdot/h;
+    g(problem.nstates*(i-1)+1+(1:problem.nstates),i+problem.nstates+1) = dfdu;
+end
+end

Easy/easydyn1.m

@@ -0,0 +1,6 @@
+function [f, dfdx, dfdxdot, dfdu] = easydyn1(x, xdot)
+
+f = x+xdot;
+dfdx = speye(length(x));
+dfdxdot = speye(length(x));
+dfdu = spalloc(length(x),1,0);

Easy/easyobj.m

@@ -0,0 +1,13 @@
+function [F, g] = easyobj(x)
+
+% f1 = sum(x.^2);
+% g(:,1) = 2*transpose(x);
+% 
+% [f, dfdx, dfdxdot, dfdu] = easydyn1(x);
+% 
+% F = [f1;transpose(f);diag(dfdx)];
+% g = [transpose(g);[dfdx,dfdxdot]]
+
+global problem
+
+[F,g] = problem.easyobj(x);

GetPareto_DW.m

@@ -0,0 +1,72 @@
+% script for running the optim.m program
+% clear all
+
+% define which movement we are working on
+movement = 'Winter/Winter_normal';
+
+% general settings
+problem.Solver = 'IPOPT';
+problem.MaxIterations = 5000;
+problem.ConstraintTol = .01;
+problem.Tol = .0001;
+problem.symmetry = 1;
+problem.discretization = 'BE';		% start with backward Euler
+problem.checkderivatives = 0;		% set this to 1 to make optim.m check derivatives
+problem.debug = 0;                  % only used in debug mode
+problem.N = 60;						% start with a coarse mesh
+problem.Printinterval = 10;
+
+% define an able-bodied model and the target gait data for the simulation
+ablemodel.parameterfile 	= 'gait2d_par.xls';
+ablemodel.type 				= 'able';
+ablemodel.datafile			= movement;
+ablemodel.Wtrack 	= 1;	% weight of tracking term in optimization objective
+ablemodel.Weffort 	= 20;	% weight of muscle effort term in optimization objective
+ablemodel.effort.fatigue 	= 0;
+ablemodel.effort.Fmaxweighted = 0;
+ablemodel.effort.exponent	= 3;
+ablemodel.Wvalve 	= 0.001;	% weight of valve operating cost in optimization objective
+ablemodel.discretization = 'euler';
+ablemodel.reducedW = 0;
+ablemodel.kneeconstraint = 1;
+ablemodel.hipconstraint = 1;
+
+% define a below knee amputee model, based on the able-bodied model
+bkamodel		= ablemodel;
+bkamodel.type	= 'bka';
+% 450 Nm/rad foot stiffness in able bodied subjects (Hansen et al., J Biomech 37: 1467–1474, 2004)
+bkamodel.anklestiffness = 600;		% stiffness (Nm/rad) of prosthetic ankle
+bkamodel.ankledamping = 15;         % damping (Nm/rad/s) of prosthetic ankle
+
+bkamodel.Wmom = 5;
+bkamodel.Weffort = 2;
+bkamodel.Wtrack = 0.1;
+
+problem.model 			= bkamodel;
+problem.initialguess 	= 'Winter/Winter_normal_result_able';%'Winter/Pareto_DW/Result_3_mom_1_effort_12_track_0.6.mat';
+problem.resultfile 		= ['Winter/Pareto_DW/Result2_3_mom_', num2str(bkamodel.Wmom), '_effort_', num2str(bkamodel.Weffort), '_track_' num2str(bkamodel.Wtrack), '.mat'];
+optim_threeobj(problem)
+
+
+
+% bkamodel.Wtrack = 1;
+% bkamodel.Weffort = 20*bkamodel.Wtrack;
+% for j = 1:10
+%     bkamodel.Wmom = j/10;
+%     problem.model 			= bkamodel;
+%     problem.initialguess 	= 'Winter\Winter_normal_result_able.mat';
+%     problem.resultfile 		= ['Winter/Pareto_DW/Result_3_mom_', num2str(bkamodel.Wmom), '_effort_', num2str(bkamodel.Weffort), '_track_' num2str(bkamodel.Wtrack), '.mat'];
+%     optim_threeobj(problem);
+% end
+% 
+% bkamodel.Wmom = 1;
+% for j = 1:9
+%     bkamodel.Wtrack = j/10;
+%     bkamodel.Weffort = 20*bkamodel.Wtrack;
+%     problem.model 			= bkamodel;
+%     problem.initialguess 	= 'Winter\Winter_normal_result_able.mat';
+%     problem.resultfile 		= ['Winter/Pareto_DW/Result_3_mom_', num2str(bkamodel.Wmom), '_effort_', num2str(bkamodel.Weffort), '_track_' num2str(bkamodel.Wtrack), '.mat'];
+%     optim_threeobj(problem);
+% end
+
+